rm(list=ls()) # clears memory 
# setwd("[Insert file path here]")
dat <- read.csv("2d_CardiacData_GoatCrossModalRecognHumans.csv", stringsAsFactors = F, 
                header = T)# load data file
head(dat,10)# view the first 10 observations


## To install and load relevant packages
# install.packages("lme4", dep=T)
# install.packages("MuMIn", dep=T)
# install.packages("lmerTest", dep=T)
# install.packages("outliers", dep=T)
# install.packages("emmeans", dep=T)
library(lme4)
library(MuMIn)
library(lmerTest)
library(outliers)
library(emmeans)


## Setting factors and checking their levels
dat$goat.id <- as.factor(dat$goat.id)
levels(dat$goat.id)

dat$goat.sex <- as.factor(dat$goat.sex)
levels(dat$goat.sex)

dat$trial <- ordered(dat$trial)
levels(dat$trial)

dat$year <- as.factor(dat$year)
levels(dat$year)

dat$congruency <- as.factor(dat$congruency)
levels(dat$congruency)

dat$human.gender <- as.factor(dat$human.gender)
levels(dat$human.gender)

dat$photo <- as.factor(dat$photo)
levels(dat$photo)

dat$responsi.photo <- as.factor(dat$responsi.photo)
levels(dat$responsi.photo)

dat$voice <- as.factor(dat$voice)
levels(dat$voice)

dat$responsi.voice <- as.factor(dat$responsi.voice)
levels(dat$responsi.voice)

dat$pb.no <- as.factor(dat$pb.no)
levels(dat$pb.no)

str(dat)# to check factors have been coded correctly

dat$obs <- 1:dim(dat)[1]


# Excluding one subject as he did not look for one congruent and one incongruent trial
dat2 <- subset(dat,goat.id!="Bill")

# Excluding observations where goats were looking before the initiation of a voice 
# playback
dat3 <- subset(dat2,obs!=57 & obs!=61 & obs!=62 & obs!=79 & obs!=78 & obs!= 85 & 
               obs!=86 & obs!=88 & obs!=93 & obs!=15 & obs!=17 & obs!=42 & obs!=48)


#############################################################################################
### HEART RATE ANALYSIS

## Full model
hr.m1 <- lmer(d.hr ~ congruency*pb.no + year + goat.sex + human.gender + responsi.photo +
              responsi.voice + prelim.dur + measure.per + (1|goat.id/trial) + (1|photo) +
              (1|voice), data=dat3)
summary(hr.m1)


## Pairwise comparisons for the interaction between congruency and playback number
emmeans(hr.m1, pairwise ~ congruency*pb.no) 


## Basic model
basic.hr.m1 <- lmer(d.hr~congruency*pb.no+(1|goat.id/trial),data=dat3)
summary(basic.hr.m1)


## Reinstating observations where goats were already looking before the start of playbacks
## to evaluate the robustness of findings
hr.m2 <- lmer(d.hr ~ congruency*pb.no + year + goat.sex + human.gender + responsi.photo +
              responsi.voice + prelim.dur + measure.per + (1|goat.id/trial) + (1|photo) +
              (1|voice), data=dat2)
summary(hr.m2)


#########################################################################################################################################
### HEART RATE VARIABILITY ANALYSIS

## Full model
hrv.m1 <- lmer(d.hrv ~ congruency*pb.no + year + goat.sex + human.gender + responsi.photo +
               responsi.voice + prelim.dur + measure.per + (1|goat.id/trial) + (1|voice) + 
               (1|photo), data=dat3)

# Removing model outliers
dat3$zscore <- scores(dat3$d.hrv, type="z", prob=0.95)
dat4 <- subset(dat3, zscore==FALSE)

# Refitting the model with the outliers removed
hrv.m2 <- lmer(d.hrv ~ congruency*pb.no + year + goat.sex + human.gender + responsi.photo +
               responsi.voice + prelim.dur + measure.per + (1|goat.id/trial) + (1|voice) +
               (1|photo), data=dat4)
summary(hrv.m2)


## Pairwise comparisons for the interaction between congruency and playback number
emmeans(hrv.m2, pairwise ~ congruency*pb.no)


## Basic model
basic.hrv.m2 <- lmer(d.hrv~congruency*pb.no+(1|goat.id/trial),data=dat4)

summary(basic.hrv.m2)


## Reinstating observations where goats were already looking before the start of playbacks
## to evaluate the robustness of findings
hrv.m3 <- lmer(d.hrv~congruency*pb.no+year+goat.sex+human.gender+responsi.photo+
               responsi.voice+prelim.dur+measure.per+(1|goat.id/trial)+(1|voice)+(1|photo),
               data=dat2)

# Removing model outliers
dat2$zscore <- scores(dat2$d.hrv, type="z", prob=0.95)
dat5 <- subset(dat2, zscore==FALSE)

hrv.m4 <- lmer(d.hrv ~ congruency*pb.no + year + goat.sex + human.gender + responsi.photo +
               responsi.voice + prelim.dur + measure.per + (1|goat.id/trial) + (1|voice) +
               (1|photo), data=dat5)
summary(hrv.m4)


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