######################################
#Load packages
install.packages("carData")
install.packages("ggplot2")                 
install.packages("ggpubr")
install.packages("car")
install.packages("lattice")
install.packages("latticeExtra")
library(lattice)
library(ggplot2)
library(tidyverse)
library(lme4)
library(car)
library(carData)
library(ggpubr)
library(latticeExtra)


##################### Wilcoxon test ######################

wilcox.test(Reversal$Pre.L.eye,
            Reversal$Post.L.eye,
            paired = TRUE,
            conf.level = 0.95)

wilcox.test(Reversal$Pre.R.eye,
            Reversal$Post.R.eye,
            paired = TRUE,
            conf.level = 0.95)

wilcox.test(Reversal$Baseline.BR,
            Reversal$Blink.rate,
            paired = TRUE,
            conf.level = 0.95)

wilcox.test(Reversal$Baseline.RMSSD,
            Reversal$RMSSD,
            paired = TRUE,
            conf.level = 0.95)


######################################

# missing values?
colSums(is.na(Reversal))


# categorical covariates balanced?
table(Reversal$Sex)

table(Reversal$Type)

table(Reversal$Background)

# OUTLIERS

# vector of continuous covariates
Var <- c("Left.eye.change", "Right.eye.change", "Blink.Rate", 
         "RMSSD","Attempts", "Sex", "Age", "Breed", "Type", 
         "Colour", "Background", "Baseline.BR", "Baseline.RMSSD")
Var

# multi-panel Cleveland dotplot

#Fig 1.
dotplot(as.matrix(as.matrix(Reversal[Var,])),
        groups=FALSE,
        strip = strip.custom(bg = 'white',
                             par.strip.text = list(cex = 1.2)),
        scales = list(x = list(relation = "free", draw = TRUE),
                      y = list(relation = "free", draw = FALSE)),
        col=1, cex  = 0.6, pch = 16,
        xlab = list(label = "Data range", cex = 1.2),
        ylab = list(label = "Data order", cex = 1.2))



######################################
#NUMBER OF ZEROS IN THE RESPONSE VARIABLE

sum(Reversal$taxa == 0)


######################################
#COLLINEARITY


Coll <- c("Left.eye.change","Right.eye.change",
        "RMSSD","Baseline.RMSSD","Blink.Rate",
          "Baseline.BR","Performance")

#Fig 2.
panel.cor <- function(x, y, digits=1, prefix="", cex.cor = 6)
{usr <- par("usr"); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r1=cor(x,y,use="pairwise.complete.obs")
r <- abs(cor(x, y,use="pairwise.complete.obs"))
txt <- format(c(r1, 0.1), digits=digits)[1]
txt <- paste(prefix, txt, sep="")
if(missing(cex.cor)) { cex <- 0.6/strwidth(txt) } else {
  cex = cex.cor}
text(0.5, 0.5, txt, cex = cex * r)}
pairs(Reversal[,Coll], lower.panel = panel.cor, cex.labels = 1.5)


#VIF
vif(glm(Performance ~  Left.eye.change + Right.eye.change +
          Blink.Rate + RMSSD + Baseline.BR +
          Baseline.RMSSD,
        family = poisson,
        data = Reversal))

######################################
#PLOT RELATIONSHIPS

# Plot response variable against  covariates

# Fig. 3.
par(mfrow=c(1,1), mar=c(5,5,1,1))
plot(Performance ~ Left.eye.change,  data = Reversal, 
     xlab = "Temperature change (Celcius)", ylab = "Performance Index",
     cex.lab = 1.8,
     pch = 16, cex = 1.3)
abline(lm(Reversal$Performance ~ Reversal$Left.eye.change),
            col = "black")

plot(Performance ~ Baseline.RMSSD,  data = Reversal, 
     xlab = "Baseline HRV (RMSSD)", ylab = "Performance Index",
     cex.lab = 1.8,
     pch = 16, cex = 1.3)
abline(lm(Reversal$Performance ~ Reversal$Baseline.RMSSD),
       col = "black")
fit1 <- lm(Performance~Baseline.RMSSD, data=Reversal)
fit2 <- lm(Performance~poly(Baseline.RMSSD,2,raw=TRUE), data=Reversal)
fit3 <- lm(Performance~poly(Baseline.RMSSD,3,raw=TRUE), data=Reversal)
fit4 <- lm(Performance~poly(Baseline.RMSSD,4,raw=TRUE), data=Reversal)
fit5 <- lm(Performance~poly(Baseline.RMSSD,5,raw=TRUE), data=Reversal)
fit6 <- lm(Performance~poly(Baseline.RMSSD,6,raw=TRUE), data=Reversal)
fit7 <- lm(Performance~poly(Baseline.RMSSD,7,raw=TRUE), data=Reversal)
fit8 <- lm(Performance~poly(Baseline.RMSSD,8,raw=TRUE), data=Reversal)

range(Reversal$Baseline.RMSSD)
x_axis <- seq(26.17,117.12,length=20)

lines(x_axis, predict(fit1, data.frame(Baseline.RMSSD=x_axis)), col='green')
lines(x_axis, predict(fit2, data.frame(Baseline.RMSSD=x_axis)), col='red')
lines(x_axis, predict(fit3, data.frame(Baseline.RMSSD=x_axis)), col='purple')
lines(x_axis, predict(fit4, data.frame(Baseline.RMSSD=x_axis)), col='blue')
lines(x_axis, predict(fit5, data.frame(Baseline.RMSSD=x_axis)), col='orange')
lines(x_axis, predict(fit6, data.frame(Baseline.RMSSD=x_axis)), col='pink')
lines(x_axis, predict(fit7, data.frame(Baseline.RMSSD=x_axis)), col='black')
lines(x_axis, predict(fit8, data.frame(Baseline.RMSSD=x_axis)), col='yellow')

summary(fit1)$adj.r.squared
summary(fit2)$adj.r.squared
summary(fit3)$adj.r.squared
summary(fit4)$adj.r.squared
summary(fit5)$adj.r.squared
summary(fit6)$adj.r.squared
summary(fit7)$adj.r.squared
summary(fit8)$adj.r.squared

###Line 6 has the best fit###

plot(Performance ~ Baseline.RMSSD,  data = Reversal, 
     xlab = "Baseline HRV (RMSSD)", ylab = "Performance Index",
     cex.lab = 1.8,
     pch = 16, cex = 1.3)
fit6 <- lm(Performance~poly(Baseline.RMSSD,6,raw=TRUE), data=Reversal)
lines(x_axis, predict(fit6, data.frame(Baseline.RMSSD=x_axis)), col='navy')




plot(Performance ~ Blink.Rate,  data = Reversal, 
     xlab = "Blink rate", ylab = "Performance index",
     cex.lab = 1.8,
     pch = 16, cex = 1.3)
abline(lm(Reversal$Performance ~ Reversal$Blink.Rate),
       col = "black")

plot(Performance ~ Baseline.BR,  data = Reversal, 
     xlab = "Baseline blink rate", ylab = "Performance index",
     cex.lab = 1.8,
     pch = 16, cex = 1.3)
abline(lm(Reversal$Performance ~ Reversal$Baseline.BR),
       col = "black")



##############################################

#Removing NAs
Reversal <- Reversal[complete.cases(Reversal), ]

dim(Reversal)

#Poisson GLM including all covariates
Pois1 <- glm(Performance ~  Left.eye.change + 
               Right.eye.change + Blink.Rate + 
               RMSSD + Baseline.RMSSD 
             + Baseline.BR,
             data = Reversal, 
             family = poisson(link = log))

summary(Pois1)

# MODEL VALIDATION

# Assess overdispersion
ods1 <- Pois1$deviance / Pois1$df.residual
ods1

#Biological covariance by adding Heart rate data- does this solve overdispersion?
Pois2 <- glm(Performance ~  Left.eye.change * 
               Right.eye.change + Blink.Rate + 
               Heart.Rate * RMSSD,
             data = Reversal, 
             family = poisson(link = log))

ods2 <- Pois2$deviance / Pois2$df.residual
ods2

#Still Overdispersed
#Use negative binomial glm instead.


#2. Zero inflation

sum(Reversal$Performance == 0)
#no zeros


#3. Influential outliers


# Fig 4. Cook's distance
par(mfrow=c(1,1), mar=c(5,5,2,2))
plot(cooks.distance(Pois1),
     xlab = "Observation", 
     ylab = "Cook's distance",
     type = "h", 
     ylim = c(0, 1.2),
     cex.lab =  1.5)
abline(h = 1, lty = 2)


#4. Non-independence

# Plot performance against horse

#Fig. 5.
par(mfrow=c(1,1), mar=c(5,5,2,2))
boxplot(Performance ~ Horse, 
        data = Reversal, 
        xlab = "Horse",
        ylab = "Performance Index",
        pch = 16, cex.lab = 1.5, col = "steelblue2")


#5. Wrong link function

# Try alternative link functions
#Identity link
Pois3 <- glm(Performance~  Left.eye.change + 
               Right.eye.change + Blink.Rate + 
               RMSSD +Baseline.RMSSD +
              Baseline.BR,
             data = Reversal, 
             family = poisson(link = identity))
summary(Pois3)

ods3 <- Pois3$deviance / Pois3$df.residual
ods3

#Square-root link
Pois4 <- glm(Performance~  Left.eye.change + 
               Right.eye.change + Blink.Rate + 
               RMSSD +Baseline.RMSSD +
               Baseline.BR,
             data = Reversal, 
             family = poisson(link = sqrt))
ods4 <- Pois4$deviance / Pois4$df.residual
ods4
summary(Pois4)



E1 <- resid(Pois1, type = "pearson")

#Fig 6.

xyplot(E1 ~ Left.eye.change, 
       data = Reversal,
       ylab = list("Pearson residuals", cex = 1.5),
       xlab = list("Left eye change", cex = 1.5),
       strip = function(bg='white', ...)
         strip.default(bg='white', ...),
       scales = list(alternating = TRUE,
                     x = list(relation = "free"),
                     y = list(relation = "same")),
       panel = function(x,y){
         panel.points(x,y, col = 1, pch = 16, cex = 1.0)
         panel.loess(x,y, col = 1, lwd = 3)})

#Fig 7. 

xyplot(E1 ~ Right.eye.change, 
       data = Reversal,
       ylab = list("Pearson residuals", cex = 1.5),
       xlab = list("Right eye change", cex = 1.5),
       strip = function(bg='white', ...)
         strip.default(bg='white', ...),
       scales = list(alternating = TRUE,
                     x = list(relation = "free"),
                     y = list(relation = "same")),
       panel = function(x,y){
         panel.points(x,y, col = 1, pch = 16, cex = 1.0)
         panel.loess(x,y, col = 1, lwd = 3)})
# Fig 8.

xyplot(E1 ~ Blink.Rate, 
       data = Reversal,
       ylab = list("Pearson residuals", cex = 1.5),
       xlab = list("Blink rate", cex = 1.5),
       strip = function(bg='white', ...)
         strip.default(bg='white', ...),
       scales = list(alternating = TRUE,
                     x = list(relation = "free"),
                     y = list(relation = "same")),
       panel = function(x,y){
         panel.points(x,y, col = 1, pch = 16, cex = 1.0)
         panel.loess(x,y, col = 1, lwd = 3)})

#Fig 9.
xyplot(E1 ~ Baseline.RMSSD, 
       data = Reversal,
       ylab = list("Pearson residuals", cex = 1.5),
       xlab = list("Baseline.RMSSD", cex = 1.5),
       strip = function(bg='white', ...)
         strip.default(bg='white', ...),
       scales = list(alternating = TRUE,
                     x = list(relation = "free"),
                     y = list(relation = "same")),
       panel = function(x,y){
         panel.points(x,y, col = 1, pch = 16, cex = 1.0)
         panel.loess(x,y, col = 1, lwd = 3)})


##negative binomial glm###

library(MASS)
nb1 <- glm.nb(Performance ~  Blink.rate + Baseline.BR +
                Left.eye.change + Right.eye.change
              +RMSSD + Baseline.RMSSD,
              data = Reversal)
summary(nb1)

nb2 <- glm.nb(Performance ~ Left.eye.change + Baseline.RMSSD,
              data = Reversal)
summary(nb2)

ods_nb <- nb1$deviance / nb1$df.residual
ods_nb



#Plot the residuals vs fitted values.
#Fig 8.7
E2 <- resid(nb2, test='pearson')
F2 <- fitted(nb2)

par(mfrow = c(2,3), mar = c(5,5,1,1))
plot(x = F2,
     y = E2,
     pch = 16, cex.lab = 1.2,
     xlab = "Fitted values",
     ylab = "Pearson residuals")
text(30, 2.1, "A", cex = 1.2)
abline(h = 0, lty = 2)
abline(v = 0, lty = 2)

#Plot residuals against all covariates in the model

plot(x = Reversal$Left.eye.change, 
     y = E2,
     xlab = "Left eye change",
     ylab = "Pearson residuals",
     pch = 16, cex.lab = 1.2,
     type = "p")
text(0.55, 2.1, "B", cex = 1.2)
abline(h = 0, lty = 2)


# Plot cook's distance to identify influential observations.


par(mfrow=c(1,1), mar=c(5,5,2,2))
plot(cooks.distance(nb1),
     xlab = "Observation", 
     ylab = "Cook's distance",
     type = "h", 
     ylim = c(0, 1.2),
     cex.lab =  1.5)
abline(h = 1, lty = 2)


# Compare fit of poisson and negative binomial models
AIC(nb2,Pois1)


# Model summary
nb1 <- glm.nb(Performance ~ Baseline.RMSSD+
                 Left.eye.change,
              data = Reversal)

summary(nb1)

nb2 <- glm.nb(Performance ~ Left.eye.change 
              + Baseline.RMSSD
              , data = Reversal)

summary(nb2)

plot(nb2)

Coll <- c("Left.eye.change","Baseline.RMSSD","Performance")

#Fig 2.
panel.cor <- function(x, y, digits=1, prefix="", cex.cor = 6)
{usr <- par("usr"); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r1=cor(x,y,use="pairwise.complete.obs")
r <- abs(cor(x, y,use="pairwise.complete.obs"))
txt <- format(c(r1, 0.1), digits=digits)[1]
txt <- paste(prefix, txt, sep="")
if(missing(cex.cor)) { cex <- 0.6/strwidth(txt) } else {
  cex = cex.cor}
text(0.5, 0.5, txt, cex = cex * r)}
pairs(Reversal[,Coll], lower.panel = panel.cor, cex.labels = 1.5)


########################################################


range(Reversal$Blink.Rate) #8.22-32.42

MyData <- expand.grid(
  Left.eye.change = mean(Reversal$Left.eye.change),
  Right.eye.change = mean(Reversal$Right.eye.change),
  Blink.Rate = mean(Reversal$Blink.Rate),
  RMSSD = mean(Reversal$RMSSD))
head(Reversal)

X <- model.matrix(Performance ~  Left.eye.change + 
                    Right.eye.change + Blink.Rate + 
                    RMSSD,
                  data = Reversal)
head(MyData)


##### Plot the variables as terciles of each other ###################

library(tidyverse)

Reversal <- read.csv("Reversal")

##### Relationship between baseline rmssd 
##### and performance, in each tercile of 
##### left eye temperature change:

Reversal %>% 
  mutate(rmssd_ntile = ntile(Left.eye.change, 3)) %>%
  ggplot(aes(x = Baseline.RMSSD, y = Performance)) + 
  geom_point() + 
  stat_smooth(se = T, method = 'lm', formula = y ~ x, fullrange = T) + 
  facet_wrap(~rmssd_ntile)


##### Relationship between left eye temp change 
##### and performance, in each tercile of 
##### baseline rmssd:

Reversal %>% 
  mutate(rmssd_ntile = ntile(Baseline.RMSSD, 3)) %>%
  ggplot(aes(x = Left.eye.change, y = Performance)) + 
  geom_point() + 
  stat_smooth(se = T, method = 'lm', formula = y ~ x, fullrange = T) + 
  facet_wrap(~rmssd_ntile)




#################################################################################

#######ANOVA#########
Reversal$group <- ordered(Reversal$group,
                         levels = c("Blink.Rate", "Baseline.BR"))

library(ggplot2)
library(ggpubr)
theme_set(theme_pubr())
hist(Reversal$Performance,
     xlab = "Performance Index",
     ylab = "Frequency",
     main = "Histogram of performance index")