#### 0. Clean up environment, download and load packages ####

rm(list = ls())

#install.packages("tidyverse")
#install.packages("readxl")
#install.packages("openxlsx")
#install.packages("performance")
#install.packages("lme4")
#install.packages("lmerTest")
#install.packages("MuMIn")
#install.packages("arm")
#install.packages("merTools")
#install.packages("cowplot")

library(tidyverse)
library(readxl)
library(openxlsx)
library(performance)
library(lme4)
library(lmerTest)
library(MuMIn)
library(arm)
library(merTools)
library(cowplot)


#### 1. Read and format data for analysis ####

# I have kept the datasets for the various indicators separate as the time points are quite different

# Read demographics info to merge with other data
dem <- read_excel("Demographics.xlsx")

#### IRT data

# IRT
IRT <- read_excel("Pain & Personality MASTER.xlsx", 
                  sheet = "IRT") %>% 
  dplyr::select(Subject, Name, IRT.BLAv, CastIRT.Av, IRT1.Av, IRT2.Av, IRT3.Av, IRT4.Av, IRT5.Av, IRT6.Av, IRT7.Av, IRT8.Av, IRT24.Av, IRT48.Av, IRT72.Av) %>% 
  rename(Horse_ID = Subject,
         IRTBaseline = IRT.BLAv,
         IRTCastration = CastIRT.Av) %>%
  left_join(., dem, by = c("Horse_ID", "Name")) %>%
  pivot_longer(., cols = contains("IRT"), names_to = "TimePoint", values_to = "IRT") %>%
  mutate(TimePoint = str_remove(TimePoint,"IRT"),
         TimePoint = str_remove(TimePoint, ".Av"))

# IRT response
IRT_resp <- read_excel("Pain & Personality MASTER.xlsx", 
                       sheet = "IRT") %>% 
  dplyr::select(Subject, IRT.BLAv, CastIRT.Av, IRT1.Av, IRT2.Av, IRT3.Av, IRT4.Av, IRT5.Av, IRT6.Av, IRT7.Av, IRT8.Av, IRT24.Av, IRT48.Av, IRT72.Av) %>%
  mutate(across(-matches('Subject'), ~. - IRT.BLAv)) %>%
  dplyr::select(-c(IRT.BLAv)) %>%
  pivot_longer(., cols = contains("IRT"), names_to = "TimePoint", values_to = "IRT_resp") %>%
  mutate(TimePoint = str_remove(TimePoint,"IRT"),
         TimePoint = str_remove(TimePoint, ".Av"),
         TimePoint= str_replace(TimePoint, "Cast", "Castration"))

IRT <- left_join(IRT, IRT_resp, by = join_by(Horse_ID == Subject, TimePoint == TimePoint))  
rm(IRT_resp)

# Difference in right and left eye temperature
IRT_Diff <- read_excel("Pain & Personality MASTER.xlsx", 
                       sheet = "IRT") %>% 
  dplyr::select("Name", "PreBLIRT.Left", "PrePTIRT.Left", "CastIRT.Left", "IRT1.Left", "IRT2.Left", "IRT3.Left", "IRT4.Left", "IRT5.Left", "IRT6.Left", "IRT7.Left", "IRT8.Left", "IRT24.Left", "IRT48.Left", "IRT72.Left" ) %>%
  mutate(IRTBaseline.Left = (PreBLIRT.Left + PrePTIRT.Left) / 2) %>%
  dplyr::select(!(contains("PreBLIRT") | contains("PrePTIRT"))) %>%
  rename_with(~str_remove(., ".Left")) %>%
  pivot_longer(., cols = contains("IRT"), names_to = "TimePoint", values_to = "IRT_Left")

IRT_Right <- read_excel("Pain & Personality MASTER.xlsx", 
                        sheet = "IRT") %>% 
  dplyr::select("Name", "PreBLIRT.Right", "PrePTIRT.Right", "CastIRT.Right", "IRT1.Right", "IRT2.Right", "IRT3.Right", "IRT4.Right", "IRT5.Right", "IRT6.Right", "IRT7.Right", "IRT8.Right", "IRT24.Right", "IRT48.Right", "IRT72.Right" ) %>%
  mutate(IRTBaseline.Right = (PreBLIRT.Right + PrePTIRT.Right) / 2) %>%
  dplyr::select(!(contains("PreBLIRT") | contains("PrePTIRT"))) %>%
  rename_with(~str_remove(., ".Right")) %>%
  pivot_longer(., cols = contains("IRT"), names_to = "TimePoint", values_to = "IRT_Right")

IRT_Diff <- left_join(IRT_Diff, IRT_Right, by = c("Name", "TimePoint")) %>%
  mutate(IRT_Diff = IRT_Right - IRT_Left) %>%
  mutate(TimePoint = str_remove(TimePoint,"IRT"),
         TimePoint = str_remove(TimePoint, ".Av"),
         TimePoint = str_replace(TimePoint, "Cast", "Castration")) %>% 
  group_by(Name) %>% 
  mutate(IRT_Diff_resp = IRT_Diff - IRT_Diff[TimePoint == "Baseline"]) %>%
  ungroup()


IRT <- left_join(IRT, dplyr::select(IRT_Diff, c("Name", "TimePoint", "IRT_Diff", "IRT_Diff_resp")), by = join_by(Name == Name, TimePoint == TimePoint))
rm(IRT_Diff, IRT_Right)

IRT$TimePoint <- factor(IRT$TimePoint, levels = unique(IRT$TimePoint))


#### HGS data
HGS <- read_excel("Pain & Personality MASTER.xlsx", 
                  sheet = "HGS") %>% 
  dplyr::select(Subject, Name, BLHGS.Av, PostCast.HGS, Post1.HGS, Post2.HGS, Post3.HGS, Post4.HGS, Post8.HGS, Post24.HGS, Post48.HGS, Post72.HGS) %>% 
  rename(Horse_ID = Subject,
         Baseline.HGS = BLHGS.Av,
         Castration.HGS = PostCast.HGS) %>%
  left_join(., dem, by = c("Horse_ID", "Name")) %>%
  pivot_longer(., cols = contains("HGS"), names_to = "TimePoint", values_to = "HGS") %>%
  mutate(TimePoint = str_remove(TimePoint,".HGS"),
         TimePoint = str_remove(TimePoint, "Post"))

HGS_resp <- read_excel("Pain & Personality MASTER.xlsx", 
                       sheet = "HGS") %>% 
  dplyr::select(Subject, BLHGS.Av, PostCast.HGS, Post1.HGS, Post2.HGS, Post3.HGS, Post4.HGS, Post8.HGS, Post24.HGS, Post48.HGS, Post72.HGS) %>%
  mutate(across(-matches('Subject'), ~. - BLHGS.Av)) %>%
  dplyr::select(-c(BLHGS.Av)) %>%
  pivot_longer(., cols = contains("HGS"), names_to = "TimePoint", values_to = "HGS_resp") %>%
  mutate(TimePoint = str_remove(TimePoint,".HGS"),
         TimePoint = str_remove(TimePoint, "Post"),
         TimePoint= str_replace(TimePoint, "Cast", "Castration"))

HGS <- left_join(HGS, HGS_resp, by = join_by(Horse_ID == Subject, TimePoint == TimePoint))  
rm(HGS_resp)

HGS$TimePoint <- factor(HGS$TimePoint, levels = unique(HGS$TimePoint))


# Maintenance and Pain Behaviour #
BEH <- read_excel("Pain & Personality MASTER.xlsx", 
                  sheet = "Etho_Dur") %>%
  dplyr::select(Name, Maintenance_PRE, Maintenance_4h, Maintenance_8h, Maintenance_24h, Maintenance_48h, Maintenance_72h) %>% 
  rename(Maintenance_Baseline = Maintenance_PRE) %>%
  left_join(., dem, by = c("Name")) %>%
  pivot_longer(., cols = contains("Maintenance"), names_to = "TimePoint", values_to = "M_BEH") %>%
  mutate(TimePoint = str_remove(TimePoint,"Maintenance_"),
         TimePoint = str_remove(TimePoint, "h"))

M_BEH_resp <- read_excel("Pain & Personality MASTER.xlsx", 
                         sheet = "Etho_Dur") %>%
  dplyr::select(Name, Maintenance_PRE, Maintenance_4h, Maintenance_8h, Maintenance_24h, Maintenance_48h, Maintenance_72h) %>% 
  mutate(across(-matches('Name'), ~. - Maintenance_PRE)) %>%
  dplyr::select(-c(Maintenance_PRE)) %>%
  pivot_longer(., cols = contains("Maintenance"), names_to = "TimePoint", values_to = "M_BEH_resp") %>%
  mutate(TimePoint = str_remove(TimePoint,"Maintenance_"),
         TimePoint = str_remove(TimePoint, "h"))

BEH <- left_join(BEH, M_BEH_resp, by = join_by(Name == Name, TimePoint == TimePoint))  
rm(M_BEH_resp)

P_BEH <- read_excel("Pain & Personality MASTER.xlsx", 
                    sheet = "Etho_Dur") %>%
  dplyr::select(Name, Pain_rel_PRE, Pain_rel_4h, Pain_rel_8h, Pain_rel_24h, Pain_rel_48h, Pain_rel_72h) %>% 
  rename(Pain_rel_Baseline = Pain_rel_PRE) %>%
  pivot_longer(., cols = contains("Pain_rel_"), names_to = "TimePoint", values_to = "P_BEH") %>%
  mutate(TimePoint = str_remove(TimePoint,"Pain_rel_"),
         TimePoint = str_remove(TimePoint, "h"))

BEH <- left_join(BEH, P_BEH, by = join_by(Name == Name, TimePoint == TimePoint))  
rm(P_BEH)

P_BEH_resp <- read_excel("Pain & Personality MASTER.xlsx", 
                         sheet = "Etho_Dur") %>%
  dplyr::select(Name, Pain_rel_PRE, Pain_rel_4h, Pain_rel_8h, Pain_rel_24h, Pain_rel_48h, Pain_rel_72h) %>% 
  mutate(across(-matches('Name'), ~. - Pain_rel_PRE)) %>%
  dplyr::select(-c(Pain_rel_PRE)) %>%
  pivot_longer(., cols = contains("Pain_rel_"), names_to = "TimePoint", values_to = "P_BEH_resp") %>%
  mutate(TimePoint = str_remove(TimePoint,"Pain_rel_"),
         TimePoint = str_remove(TimePoint, "h"))

BEH <- left_join(BEH, P_BEH_resp, by = join_by(Name == Name, TimePoint == TimePoint))  
rm(P_BEH_resp)

BEH$TimePoint <- factor(BEH$TimePoint, levels = unique(BEH$TimePoint))


# Cortisol #
CORT <- read_excel("Pain & Personality MASTER.xlsx", 
                   sheet = "Cortisol") %>%
  dplyr::select(Name, `BLAv. Cort`, PreCast.Cortisol, PostCast.Cortisol, Post2.Cortisol, Post3.Cortisol, `Post4. Cortisol`, `Post8. Cortisol`, `Post24. Cortisol`, `Post48. Cortisol`, `Post72. Cortisol`) %>% 
  rename_with(~str_remove(., " ")) %>%
  rename(Baseline.Cortisol = BLAv.Cort) %>%
  left_join(., dem, by = c("Name")) %>%
  pivot_longer(., cols = contains("Cortisol"), names_to = "TimePoint", values_to = "CORT") %>%
  mutate(TimePoint = str_remove(TimePoint,".Cortisol"),
         TimePoint = str_remove(TimePoint, "Post"),
         TimePoint= str_replace(TimePoint, "Cast", "Castration"))

CORT_resp <- read_excel("Pain & Personality MASTER.xlsx", 
                        sheet = "Cortisol") %>%
  dplyr::select(Name, `BLAv. Cort`, PreCast.Cortisol, PostCast.Cortisol, Post2.Cortisol, Post3.Cortisol, `Post4. Cortisol`, `Post8. Cortisol`, `Post24. Cortisol`, `Post48. Cortisol`, `Post72. Cortisol`) %>% 
  rename_with(~str_remove(., " ")) %>%
  mutate(across(-matches('Name'), ~. - BLAv.Cort)) %>%
  dplyr::select(-c(BLAv.Cort)) %>%
  pivot_longer(., cols = contains("Cortisol"), names_to = "TimePoint", values_to = "CORT_resp") %>%
  mutate(TimePoint = str_remove(TimePoint,".Cortisol"),
         TimePoint = str_remove(TimePoint, "Post"),
         TimePoint= str_replace(TimePoint, "Cast", "Castration"))

CORT <- left_join(CORT, CORT_resp, by = join_by(Name == Name, TimePoint == TimePoint)) 
rm(CORT_resp)

CORT$TimePoint <- factor(CORT$TimePoint, levels = unique(CORT$TimePoint))

rm(dem)


# All variables at TimePoint4 #

ALL_T4 <- Reduce(function(dtf1,dtf2) left_join(dtf1,dtf2,by="Horse_ID"), 
                 list(
                   filter(IRT, TimePoint == "4"),
                   filter(CORT, TimePoint == "4") %>% dplyr::select("Horse_ID", "CORT", "CORT_resp"),
                   filter(HGS, TimePoint == "4") %>% dplyr::select("Horse_ID", "HGS", "HGS_resp"),
                   filter(BEH, TimePoint == "4") %>% dplyr::select("Horse_ID", "M_BEH", "M_BEH_resp", "P_BEH", "P_BEH_resp")))



#### 2. Data analysis ####


#### 2.1. AIM 1: Is there a significant impact of time point on the ABSOLUTE VALUE of the physiological and behavioural variables? ####

# 2.1.1. Impact of timepoints on IRT 
IRT_mod <- lmer(data = IRT, IRT ~ TimePoint + Breed + (1|Horse_ID) + (1|Farm) + (1|Veterinarian))
check_model(IRT_mod)
summary(IRT_mod)
r.squaredGLMM(IRT_mod)

# IRT is significantly different from baseline from immediately post castration to 5hrs post castration
# Restrict repeatability of IRT response to those time points


# 2.1.2. Impact of timepoint on IRT Discrepancy
IRT_Diff_mod <- lmer(data = IRT, IRT_Diff ~ TimePoint + Breed + (1|Horse_ID) + (1|Farm) + (1|Veterinarian))
check_model(IRT_mod)
summary(IRT_Diff_mod)
r.squaredGLMM(IRT_Diff_mod)

# IRT Discrepancy is actually not significantly different from baseline at any point
# Use all time points in repeatability analysis of IRT Discrepancy response


# 2.1.3. Impact of time point on HGS 
HGS_mod <- lmer(data = HGS, HGS ~ TimePoint + Breed + (1|Horse_ID) + (1|Farm) + (1|Veterinarian))
check_model(HGS_mod)
summary(HGS_mod)
r.squaredGLMM(HGS_mod)

# HGS is significantly different from baseline from 1hr post-castration to 48hr post castration
# Restrict repeatability analysis of HGS response to those time points


# 2.1.4. Impact of time point on maintenance behaviour
MBEH_mod <- lmer(data = BEH, M_BEH ~ TimePoint + Breed + (1|Horse_ID) + (1|Farm) + (1|Veterinarian))
check_model(MBEH_mod)
plot(MBEH_mod) # The residual plot looks a little wonky but I think it's mostly just a floor effect due to the hard boudary at 0
summary(MBEH_mod)
r.squaredGLMM(MBEH_mod)

# Maintenance behaviour is significantly different from baseline at all time points
# Use all time points in repeatability analysis of maintenance behaviour response


# 2.1.5. Impact of time point on pain behaviour
# (i.e. is expression of pain behaviour significantly different from baseline at the various points post-castration)
PBEH_mod <- lmer(data = BEH, P_BEH ~ TimePoint + Breed + (1|Horse_ID) + (1|Farm) + (1|Veterinarian))
check_model(PBEH_mod)
plot(PBEH_mod) # The residual plot looks very wonky but I think it's mostly just a floor effect due to the hard boudary at 0 (of which there are many)
summary(PBEH_mod)
r.squaredGLMM(PBEH_mod)

# Pain behaviour is actually not significantly different from baseline at any point
# Use all time points in repeatability analysis of pain behaviour response


# 2.1.6. Impact of time point on cortisol
# (i.e. is expression of cortisol significantly different from baseline at the various points post-castration)
CORT_mod <- lmer(data = CORT, CORT ~ TimePoint + Breed + (1|Horse_ID) + (1|Farm) + (1|Veterinarian))
check_model(CORT_mod)
plot(CORT_mod) 
# What marvellous heteroskedasticity
# Looking at the graph, timepoints where cortisol is highest are also those where it's a lot more variable, so I'm guessing that's where this comes from
summary(CORT_mod)
r.squaredGLMM(CORT_mod)

# Cortisol significantly different from baseline at 2, 3, 4, 8 and 24hrs post castration
# Restrict repeatability analysis of cortisol response to these points




#### AIM 3: Does personality (Neuroticism and Extroversion) impact the RESPONSE to castration (TimePoint - Baseline) of the physiological and behavioural variables? ####


# 2.2.1. Impact of personality on IRT response

IRT_resp_mod <- lmer(data = IRT[which(IRT$TimePoint %in% c("Castration", "1", "2", "3", "4", "5")),], IRT_resp ~ TimePoint + Breed + scale(Neuroticism) + scale(Extroversion) + (1|Horse_ID) + (1|Farm) + (1|Veterinarian))
check_model(IRT_resp_mod)
summary(IRT_resp_mod)
r.squaredGLMM(IRT_resp_mod)


# 2.2.2. Impact of personality on IRT Discrepancy RESPONSE
IRT_Diff_resp_mod <- lmer(data = IRT, IRT_Diff_resp ~ TimePoint + Breed + scale(Neuroticism) + scale(Extroversion) + (1|Horse_ID) + (1|Farm) + (1|Veterinarian))
check_model(IRT_Diff_resp_mod)
summary(IRT_Diff_resp_mod)
r.squaredGLMM(IRT_Diff_resp_mod)


# 2.2.3. Impact of personality on HGS RESPONSE
HGS_resp_mod <- lmer(data = HGS[which(HGS$TimePoint %in% c("Castration", "1", "2", "3", "4", "8", "24", "48")),], HGS_resp ~ TimePoint + Breed + scale(Neuroticism) + scale(Extroversion) + (1|Horse_ID) + (1|Farm) + (1|Veterinarian))
check_model(HGS_resp_mod)
summary(HGS_resp_mod)
r.squaredGLMM(HGS_resp_mod)


# 2.2.4. Impact of personality on Maintenance Behaviour RESPONSE
MBEH_resp_mod <- lmer(data = BEH, M_BEH_resp ~ TimePoint + Breed + scale(Neuroticism) + scale(Extroversion) + (1|Horse_ID) + (1|Farm) + (1|Veterinarian))
check_model(MBEH_resp_mod)
summary(MBEH_resp_mod)
r.squaredGLMM(MBEH_resp_mod)


# 2.2.5. Impact of personality on Pain Behaviour RESPONSE
PBEH_resp_mod <- lmer(data = BEH, P_BEH_resp ~ TimePoint + Breed + scale(Neuroticism) + scale(Extroversion) + (1|Horse_ID) + (1|Farm) + (1|Veterinarian))
check_model(PBEH_resp_mod)
summary(PBEH_resp_mod)
r.squaredGLMM(PBEH_resp_mod)


# 2.2.6. Impact of personality on Cortisol RESPONSE
CORT_resp_mod <- lmer(data = CORT[which(CORT$TimePoint %in% c("2", "3", "4", "8", "24")),], CORT_resp ~ TimePoint + Breed + scale(Neuroticism) + scale(Extroversion) + (1|Horse_ID) + (1|Farm) + (1|Veterinarian))
check_model(CORT_resp_mod)
summary(CORT_resp_mod)
r.squaredGLMM(CORT_resp_mod)



#### AIM 2: Quantifying inter-individual variation in response to castration + variation due to Farm and Vet ####

# 2.3.1. Repeatability of IRT RESPONSE (%age of remaining variance explained by Horse_ID, Farm and Vet)
IRT_resp_Rep_Horse <- VarCorr(IRT_resp_mod)$"Horse_ID"[1] /
  (VarCorr(IRT_resp_mod)$"Horse_ID"[1] +
     VarCorr(IRT_resp_mod)$"Farm"[1] +
     VarCorr(IRT_resp_mod)$"Veterinarian"[1] +
     attr(VarCorr(IRT_resp_mod), "sc")^2)
IRT_resp_Rep_Farm <- VarCorr(IRT_resp_mod)$"Farm"[1] /
  (VarCorr(IRT_resp_mod)$"Horse_ID"[1] +
     VarCorr(IRT_resp_mod)$"Farm"[1] +
     VarCorr(IRT_resp_mod)$"Veterinarian"[1] +
     attr(VarCorr(IRT_resp_mod), "sc")^2)
IRT_resp_Rep_Vet <- VarCorr(IRT_resp_mod)$"Veterinarian"[1] /
  (VarCorr(IRT_resp_mod)$"Horse_ID"[1] +
     VarCorr(IRT_resp_mod)$"Farm"[1] +
     VarCorr(IRT_resp_mod)$"Veterinarian"[1] +
     attr(VarCorr(IRT_resp_mod), "sc")^2)

RESPONSE_REPS <- t(c(IRT_resp_Rep_Horse, IRT_resp_Rep_Farm, IRT_resp_Rep_Vet))
rm(IRT_resp_Rep_Horse, IRT_resp_Rep_Farm, IRT_resp_Rep_Vet)


# 2.3.2. Repeatability of IRT Discrepancy RESPONSE (%age of remaining variance explained by Horse_ID, Farm and Vet)
IRT_Diff_resp_Rep_Horse <- VarCorr(IRT_Diff_resp_mod)$"Horse_ID"[1] /
  (VarCorr(IRT_Diff_resp_mod)$"Horse_ID"[1] +
     VarCorr(IRT_Diff_resp_mod)$"Farm"[1] +
     VarCorr(IRT_Diff_resp_mod)$"Veterinarian"[1] +
     attr(VarCorr(IRT_Diff_resp_mod), "sc")^2)
IRT_Diff_resp_Rep_Farm <- VarCorr(IRT_Diff_resp_mod)$"Farm"[1] /
  (VarCorr(IRT_Diff_resp_mod)$"Horse_ID"[1] +
     VarCorr(IRT_Diff_resp_mod)$"Farm"[1] +
     VarCorr(IRT_Diff_resp_mod)$"Veterinarian"[1] +
     attr(VarCorr(IRT_Diff_resp_mod), "sc")^2)
IRT_Diff_resp_Rep_Vet <- VarCorr(IRT_Diff_resp_mod)$"Veterinarian"[1] /
  (VarCorr(IRT_Diff_resp_mod)$"Horse_ID"[1] +
     VarCorr(IRT_Diff_resp_mod)$"Farm"[1] +
     VarCorr(IRT_Diff_resp_mod)$"Veterinarian"[1] +
     attr(VarCorr(IRT_Diff_resp_mod), "sc")^2)


RESPONSE_REPS <- rbind(RESPONSE_REPS, t(c(IRT_Diff_resp_Rep_Horse, IRT_Diff_resp_Rep_Farm, IRT_Diff_resp_Rep_Vet)))
rm(IRT_Diff_resp_Rep_Horse, IRT_Diff_resp_Rep_Farm, IRT_Diff_resp_Rep_Vet)


# 2.3.3. Repeatability of HGS RESPONSE (%age of remaining variance explained by Horse_ID, Farm and Vet)
HGS_resp_Rep_Horse <- VarCorr(HGS_resp_mod)$"Horse_ID"[1] /
  (VarCorr(HGS_resp_mod)$"Horse_ID"[1] +
     VarCorr(HGS_resp_mod)$"Farm"[1] +
     VarCorr(HGS_resp_mod)$"Veterinarian"[1] +
     attr(VarCorr(HGS_resp_mod), "sc")^2)
HGS_resp_Rep_Farm <- VarCorr(HGS_resp_mod)$"Farm"[1] /
  (VarCorr(HGS_resp_mod)$"Horse_ID"[1] +
     VarCorr(HGS_resp_mod)$"Farm"[1] +
     VarCorr(HGS_resp_mod)$"Veterinarian"[1] +
     attr(VarCorr(HGS_resp_mod), "sc")^2)
HGS_resp_Rep_Vet <- VarCorr(HGS_resp_mod)$"Veterinarian"[1] /
  (VarCorr(HGS_resp_mod)$"Horse_ID"[1] +
     VarCorr(HGS_resp_mod)$"Farm"[1] +
     VarCorr(HGS_resp_mod)$"Veterinarian"[1] +
     attr(VarCorr(HGS_resp_mod), "sc")^2)


RESPONSE_REPS <- rbind(RESPONSE_REPS, t(c(HGS_resp_Rep_Horse, HGS_resp_Rep_Farm, HGS_resp_Rep_Vet)))
rm(HGS_resp_Rep_Horse, HGS_resp_Rep_Farm, HGS_resp_Rep_Vet)


# 2.3.4. Repeatability of Maintenance Behaviour RESPONSE (%age of remaining variance explained by Horse_ID, Farm and Vet)
MBEH_resp_Rep_Horse <- VarCorr(MBEH_resp_mod)$"Horse_ID"[1] /
  (VarCorr(MBEH_resp_mod)$"Horse_ID"[1] +
     VarCorr(MBEH_resp_mod)$"Farm"[1] +
     VarCorr(MBEH_resp_mod)$"Veterinarian"[1] +
     attr(VarCorr(MBEH_resp_mod), "sc")^2)
MBEH_resp_Rep_Farm <- VarCorr(MBEH_resp_mod)$"Farm"[1] /
  (VarCorr(MBEH_resp_mod)$"Horse_ID"[1] +
     VarCorr(MBEH_resp_mod)$"Farm"[1] +
     VarCorr(MBEH_resp_mod)$"Veterinarian"[1] +
     attr(VarCorr(MBEH_resp_mod), "sc")^2)
MBEH_resp_Rep_Vet <- VarCorr(MBEH_resp_mod)$"Veterinarian"[1] /
  (VarCorr(MBEH_resp_mod)$"Horse_ID"[1] +
     VarCorr(MBEH_resp_mod)$"Farm"[1] +
     VarCorr(MBEH_resp_mod)$"Veterinarian"[1] +
     attr(VarCorr(MBEH_resp_mod), "sc")^2)

RESPONSE_REPS <- rbind(RESPONSE_REPS, t(c(MBEH_resp_Rep_Horse, MBEH_resp_Rep_Farm, MBEH_resp_Rep_Vet)))
rm(MBEH_resp_Rep_Horse, MBEH_resp_Rep_Farm, MBEH_resp_Rep_Vet)


# 2.3.5. Repeatability of Pain Behaviour RESPONSE (%age of remaining variance explained by Horse_ID, Farm and Vet)
PBEH_resp_Rep_Horse <- VarCorr(PBEH_resp_mod)$"Horse_ID"[1] /
  (VarCorr(PBEH_resp_mod)$"Horse_ID"[1] +
     VarCorr(PBEH_resp_mod)$"Farm"[1] +
     VarCorr(PBEH_resp_mod)$"Veterinarian"[1] +
     attr(VarCorr(PBEH_resp_mod), "sc")^2)
PBEH_resp_Rep_Farm <- VarCorr(PBEH_resp_mod)$"Farm"[1] /
  (VarCorr(PBEH_resp_mod)$"Horse_ID"[1] +
     VarCorr(PBEH_resp_mod)$"Farm"[1] +
     VarCorr(PBEH_resp_mod)$"Veterinarian"[1] +
     attr(VarCorr(PBEH_resp_mod), "sc")^2)
PBEH_resp_Rep_Vet <- VarCorr(PBEH_resp_mod)$"Veterinarian"[1] /
  (VarCorr(PBEH_resp_mod)$"Horse_ID"[1] +
     VarCorr(PBEH_resp_mod)$"Farm"[1] +
     VarCorr(PBEH_resp_mod)$"Veterinarian"[1] +
     attr(VarCorr(PBEH_resp_mod), "sc")^2)

RESPONSE_REPS <- rbind(RESPONSE_REPS, t(c(PBEH_resp_Rep_Horse, PBEH_resp_Rep_Farm, PBEH_resp_Rep_Vet)))
rm(PBEH_resp_Rep_Horse, PBEH_resp_Rep_Farm, PBEH_resp_Rep_Vet)


# 2.3.6. Repeatability of Cortisol RESPONSE (%age of remaining variance explained by Horse_ID, Farm and Vet)
CORT_resp_Rep_Horse <- VarCorr(CORT_resp_mod)$"Horse_ID"[1] /
  (VarCorr(CORT_resp_mod)$"Horse_ID"[1] +
     VarCorr(CORT_resp_mod)$"Farm"[1] +
     VarCorr(CORT_resp_mod)$"Veterinarian"[1] +
     attr(VarCorr(CORT_resp_mod), "sc")^2)
CORT_resp_Rep_Farm <- VarCorr(CORT_resp_mod)$"Farm"[1] /
  (VarCorr(CORT_resp_mod)$"Horse_ID"[1] +
     VarCorr(CORT_resp_mod)$"Farm"[1] +
     VarCorr(CORT_resp_mod)$"Veterinarian"[1] +
     attr(VarCorr(CORT_resp_mod), "sc")^2)
CORT_resp_Rep_Vet <- VarCorr(CORT_resp_mod)$"Veterinarian"[1] /
  (VarCorr(CORT_resp_mod)$"Horse_ID"[1] +
     VarCorr(CORT_resp_mod)$"Farm"[1] +
     VarCorr(CORT_resp_mod)$"Veterinarian"[1] +
     attr(VarCorr(CORT_resp_mod), "sc")^2)

RESPONSE_REPS <- rbind(RESPONSE_REPS, t(c(CORT_resp_Rep_Horse, CORT_resp_Rep_Farm, CORT_resp_Rep_Vet)))
rm(CORT_resp_Rep_Horse, CORT_resp_Rep_Farm, CORT_resp_Rep_Vet)



#### AIM 4: Does HGS predict IRT and Cortisol? ####

hist(ALL_T4$HGS); shapiro.test(ALL_T4$HGS) #All good

hist(ALL_T4$IRT); shapiro.test(ALL_T4$IRT) #All good
IRT_corr <- cor.test(ALL_T4$IRT, ALL_T4$HGS)

hist(ALL_T4$CORT); shapiro.test(ALL_T4$CORT) #1 outlier and so not normal
hist(ALL_T4$CORT[which(ALL_T4$CORT <3)]); shapiro.test(ALL_T4$CORT[which(ALL_T4$CORT <3)]); shapiro.test(ALL_T4$HGS[which(ALL_T4$CORT <3)])  #removing the outlier solves the problem
CORT_corr <- cor.test(ALL_T4$CORT[which(ALL_T4$CORT <3)], ALL_T4$HGS[which(ALL_T4$CORT <3)])

# Plot the correlations

IRT_corr_plot <- ggplot(data = ALL_T4, aes(x = HGS, y = IRT)) +
  geom_point(aes(colour = Breed, shape = Farm), size = 2) +
  geom_smooth(method = lm, se = F, colour = "black") +
  theme_bw() +
  theme(legend.position = "none") +
  ylab("IRT (°C)")

CORT_corr_plot <- ggplot(data = ALL_T4, aes(x = HGS, y = CORT)) +
  geom_point(aes(colour = Breed, shape = Farm), size = 2) +
  geom_smooth(method = lm, se = F, colour = "black") +
  theme_bw() +
  ylab("Cortisol")

plot_grid(IRT_corr_plot, CORT_corr_plot, ncol = 2, rel_widths = c(1,1.4))


#### 3. Format output to export as .xlsx and save ####

#### Function to format model output of interest (Fixed effects, random effect variance/covariance table, sample sizes, R squared) ####

format_output <- function(model) {
  Fix <- as.data.frame(summary(model)$coefficients)
  colnames(Fix) <- c("Estimate", "Std_Error", "Deg_Freedom", "t_Value", "p_Value")
  Fix <- cbind(Fixed_Effect = rownames(Fix), data.frame(Fix, row.names=NULL))
  Fix$Sig = ifelse(Fix$p_Value > 0.1, "", 
                   ifelse(Fix$p_Value <= 0.1 & Fix$p_Value > 0.05, ".", 
                          ifelse(Fix$p_Value <= 0.05 & Fix$p_Value > 0.01, "*",
                                 ifelse(Fix$p_Value <= 0.01 & Fix$p_Value > 0.001, "**", "***"))))
  Fix <- Fix %>% mutate(across(where(is.numeric), \(x) round(x, 3)))
  
  Rand <- as.data.frame(summary(model)$varcor) %>% dplyr::select(grp, vcov, sdcor) %>% rename(Variable = grp, Variance = vcov, Std_Dev = sdcor) %>% mutate(across(where(is.numeric), \(x) round(x, 3)))
  
  N <- as.data.frame(summary(model)$ngrp) %>% rename(Sample_Size = 1)
  
  R_Squared <- as.data.frame(r.squaredGLMM(model)) %>% mutate(across(where(is.numeric), \(x) round(x, 3)))
  
  return(list(Fix = Fix, Rand = Rand, N = N, R_Squared = R_Squared))
}


#### 3.1. AIM 1: Format and export output to Excel ####

# 3.1.1. Format IRT output
IRT_output <- format_output(IRT_mod)
rm(IRT_mod)

# 3.1.2. IRT Discrepancy
IRT_Diff_output <- format_output(IRT_Diff_mod)
rm(IRT_Diff_mod)


# 3.1.3. HGS
HGS_output <- format_output(HGS_mod)
rm(HGS_mod)

# 3.1.4. Maintenance Behaviour
MBEH_output <- format_output(MBEH_mod)
rm(MBEH_mod)

# 3.1.5. Pain Behaviour
PBEH_output <- format_output(PBEH_mod)
rm(PBEH_mod)

# 3.1.6. Cortisol
CORT_output <- format_output(CORT_mod)
rm(CORT_mod)

# 3.1.7. Export to Excel
Aim1 = createWorkbook()
IRT_sheet = addWorksheet(Aim1, "IRT")
IRT_Diff_sheet = addWorksheet(Aim1, "IRT Discrepancy")
HGS_Sheet = addWorksheet(Aim1, "HGS")
MBEH_Sheet = addWorksheet(Aim1, "Maintenance Behaviour")
PBEH_Sheet = addWorksheet(Aim1, "Pain Behaviour")
CORT_Sheet = addWorksheet(Aim1, "Cortisol")

start.cols = c(1, 9, 12, 14)

walk2(list(IRT_output$Fix, IRT_output$Rand, IRT_output$N, IRT_output$R_Squared), start.cols, 
      ~writeData(Aim1, IRT_sheet, .x, startCol=.y))

walk2(list(IRT_Diff_output$Fix, IRT_Diff_output$Rand, IRT_Diff_output$N, IRT_Diff_output$R_Squared), start.cols, 
      ~writeData(Aim1, IRT_Diff_sheet, .x, startCol=.y))

walk2(list(HGS_output$Fix, HGS_output$Rand, HGS_output$N, HGS_output$R_Squared), start.cols, 
      ~writeData(Aim1, HGS_Sheet, .x, startCol=.y))

walk2(list(MBEH_output$Fix, MBEH_output$Rand, MBEH_output$N, MBEH_output$R_Squared), start.cols, 
      ~writeData(Aim1, MBEH_Sheet, .x, startCol=.y))

walk2(list(PBEH_output$Fix, PBEH_output$Rand, PBEH_output$N, PBEH_output$R_Squared), start.cols, 
      ~writeData(Aim1, PBEH_Sheet, .x, startCol=.y))

walk2(list(CORT_output$Fix, CORT_output$Rand, CORT_output$N, CORT_output$R_Squared), start.cols, 
      ~writeData(Aim1, CORT_Sheet, .x, startCol=.y))

saveWorkbook(Aim1, "Castration Data Analysis AIM 1 Output (Impact of timepoint on absolute value of physiological and behavioural variables).xlsx")



#### AIM 2: Format and export output to Excel ####

# 3.2.1. IRT RESPONSE
IRT_resp_output <- format_output(IRT_resp_mod)
rm(IRT_resp_mod)

# 3.2.2. IRT Discrepancy RESPONSE
IRT_Diff_resp_output <- format_output(IRT_Diff_resp_mod)
rm(IRT_Diff_resp_mod)

# 3.2.3. HGS RESPONSE
HGS_resp_output <- format_output(HGS_resp_mod)
rm(HGS_resp_mod)

# 3.2.4. Maintenance Behaviour RESPONSE
MBEH_resp_output <- format_output(MBEH_resp_mod)
rm(MBEH_resp_mod)

# 3.2.5. Pain Behaviour RESPONSE
PBEH_resp_output <- format_output(PBEH_resp_mod)
rm(PBEH_resp_mod)

# 3.2.6. Cortisol RESPONSE
CORT_resp_output <- format_output(CORT_resp_mod)
rm(CORT_resp_mod)

# 3.1.7. Export to Excel
Aim2 = createWorkbook()
IRT_sheet = addWorksheet(Aim2, "IRT")
IRT_Diff_sheet = addWorksheet(Aim2, "IRT Discrepancy")
HGS_Sheet = addWorksheet(Aim2, "HGS")
MBEH_Sheet = addWorksheet(Aim2, "Maintenance Behaviour")
PBEH_Sheet = addWorksheet(Aim2, "Pain Behaviour")
CORT_Sheet = addWorksheet(Aim2, "Cortisol")

start.cols = c(1, 9, 12, 14)

walk2(list(IRT_resp_output$Fix, IRT_resp_output$Rand, IRT_resp_output$N, IRT_resp_output$R_Squared), start.cols, 
      ~writeData(Aim2, IRT_sheet, .x, startCol=.y))

walk2(list(IRT_Diff_resp_output$Fix, IRT_Diff_resp_output$Rand, IRT_Diff_resp_output$N, IRT_Diff_resp_output$R_Squared), start.cols, 
      ~writeData(Aim2, IRT_Diff_sheet, .x, startCol=.y))

walk2(list(HGS_resp_output$Fix, HGS_resp_output$Rand, HGS_resp_output$N, HGS_resp_output$R_Squared), start.cols, 
      ~writeData(Aim2, HGS_Sheet, .x, startCol=.y))

walk2(list(MBEH_resp_output$Fix, MBEH_resp_output$Rand, MBEH_resp_output$N, MBEH_resp_output$R_Squared), start.cols, 
      ~writeData(Aim2, MBEH_Sheet, .x, startCol=.y))

walk2(list(PBEH_resp_output$Fix, PBEH_resp_output$Rand, PBEH_resp_output$N, PBEH_resp_output$R_Squared), start.cols, 
      ~writeData(Aim2, PBEH_Sheet, .x, startCol=.y))

walk2(list(CORT_resp_output$Fix, CORT_resp_output$Rand, CORT_resp_output$N, CORT_resp_output$R_Squared), start.cols, 
      ~writeData(Aim2, CORT_Sheet, .x, startCol=.y))

saveWorkbook(Aim2, "Castration Data Analysis AIM 2 Output (Impact of personality on response to castration).xlsx")




#### 3.3. AIM 3: Format and export output to Excel ####

# 3.3.1. Format
RESPONSE_REPS <- as.data.frame(RESPONSE_REPS)
colnames(RESPONSE_REPS) <- c("Horse_Adjusted_Repeatability", "Farm_Adjusted_Repeatability", "Vet_Adjusted_Repeatability")
RESPONSE_REPS <- cbind(c("IRT", "IRT Discrepancy" , "HGS", "Maintenance Behaviour", "Pain Behaviour", "Cortisol"), RESPONSE_REPS)
RESPONSE_REPS <- RESPONSE_REPS %>% mutate(across(where(is.numeric), \(x) round(x, 3)))


# 3.3.2. Export to Excel ####
Aim3 = createWorkbook()
REP_Sheet = addWorksheet(Aim3, "Repeatabilities of response")

walk2(list(RESPONSE_REPS), 1, 
      ~writeData(Aim3, REP_Sheet, .x, startCol=.y))

saveWorkbook(Aim3, "Castration Data Analysis Output AIM 3 (Repeatability of physiological and behavioural responses).xlsx")



#### 3.4. AIM 4: Format and export output to Excel ####
Aim4 = createWorkbook()
REP_Sheet = addWorksheet(Aim4, "HGS vs. Physiology correlations")

corr_output <- as.data.frame(
  rbind(
    c("IRT", round(IRT_corr$statistic, 3), round(IRT_corr$parameter, 3), round(IRT_corr$estimate, 3), round(IRT_corr$p.value, 3)), 
    c("CORT", round(CORT_corr$statistic, 3), round(CORT_corr$parameter, 3), round(CORT_corr$estimate, 3), round(CORT_corr$p.value, 3))))

colnames(corr_output) <- c("Variable", "Test Statistic", "Degrees of Freedom", "Correlations Estimate (r)", "P value")

walk2(list(corr_output), 1, 
      ~writeData(Aim4, REP_Sheet, .x, startCol=.y))

saveWorkbook(Aim4, "Castration Data Analysis Output AIM 4 (Relationship between HGS and physiological variables).xlsx")