#R v3.5.0, RStudio v0.8.3 code
#The purpose of this script is to analyse data and create figures for the manuscript 'A colonial-nesting seabird shows limited heart rate responses to natural variation in threats of polar bears" Geldart et al. submitted for publication

#Last updated August 2023 by E.A.G.

#Load packages
library(dplyr) 
require(tidyverse)
require(lme4)
require(lmerTest) 
require(ggeffects) 
require(ggplot2)
require(sjPlot)
require(car)
require(MuMIn)
require(effectsize)
require(broom)
library(lsmeans)
require(arm)

setwd("~/Documents/")
NaturalPredation <- read.csv(file.choose())

#rename parameters that will be used in analysis
NaturalPredation<-NaturalPredation%>%
  rename(eider = Recorder.ID,
         interval = Interval.Number,
         predation_HR = Total.Count,
         bear_event = Bear.Event.b,
         bear_distance = Distance.Bear.Nest..m.,
         difference_HR = HR.Differences,
         duration_viewshed_area = Duration.Viewshed.Area,
         duration_viewshed_scale = Duration.Viewshed.Area.Scaled,
         sample = Sample.Number,
         incubation_stage = incubation.stage.sample,
         temp = Temp..Celsius.,
         wind_speed = Wind.Speed..m.s.,
         intervalB = Interval.Number.B,
         duration = Duration.Bear.Viewshed..min.,
         area = Viewshed.Area..m2.)


#Double check that these variables are in the format I want
str(NaturalPredation)

#Change my random variables to factors
NaturalPredation$eider<-as.factor(NaturalPredation$eider)  
NaturalPredation$interval<-as.factor(NaturalPredation$interval)
NaturalPredation$intervalB<-as.factor(NaturalPredation$intervalB)
NaturalPredation$bear_event<-as.factor(NaturalPredation$bear_event)
NaturalPredation$sample<-as.factor(NaturalPredation$sample)

#Recheck format
str(NaturalPredation)
attach(NaturalPredation)

#Get summary statistics for each fixed effect included in analysis
summary(NaturalPredation$difference_HR)
sd(NaturalPredation$difference_HR)

summary(NaturalPredation$wind_speed)
sd(NaturalPredation$wind_speed)

summary(NaturalPredation$bear_distance)
sd(NaturalPredation$bear_distance)

summary(NaturalPredation$bear_distance)
sd(NaturalPredation$bear_distance)

summary(NaturalPredation$bear_duration)
sd(NaturalPredation$bear_duration)

summary(NaturalPredation$temp)
sd(NaturalPredation$temp)

summary(NaturalPredation$incubation_stage)
sd(NaturalPredation$incubation_stage)

summary(NaturalPredation$viewshed_area)
sd(NaturalPredation$viewshed_area)


#Run full model: 

m1FULL = lmer(difference_HR ~ bear_distance * duration_viewshed_scale + wind_speed + temp + incubation_stage + (1|eider:bear_event) + (1|eider:bear_event:sample), data = NaturalPredation, REML=TRUE) 
summary(m1FULL)

#Check model assumptions
qqnorm(residuals(m1FULL))
qqline(residuals(m1FULL))
hist(residuals(m1FULL))
boxplot(residuals(m1FULL))
shapiro.test(residuals(m1FULL))

#Data skewed, so log-transform data
NaturalPredation$difference_HR_log <- (log10(NaturalPredation$difference_HR + 10))

m1FULL = lmer(difference_HR_log ~ bear_distance * duration_viewshed_scale + wind_speed + temp + incubation_stage + (1|eider:bear_event) + (1|eider:bear_event:sample), data = NaturalPredation, REML=TRUE) 
summary(m1FULL)

#Recheck model assumptions with log transformed data
qqnorm(residuals(m1FULL))
qqline(residuals(m1FULL))
hist(residuals(m1FULL))
boxplot(residuals(m1FULL))
shapiro.test(residuals(m1FULL))

#Then run it with the random slope for distance
#Run REML=TRUE since I want to include random effect estimation to compare the 2 models that just differ in a random effect.

m1FULLslope = lmer(difference_HR_log ~ bear_distance * duration_viewshed_scale + wind_speed + temp + incubation_stage + (bear_distance|eider) + (1|eider:bear_event) + (1|eider:bear_event:sample), data = NaturalPredation, REML=TRUE) 
summary(m1FULLslope)

#Compare model with and without random slope for distance
anova(m1FULL, m1FULLslope)
#Slope does not improve fit of model

#Since slope is not significant, proceed with backward stepwise elimination removing non-significant terms, and using Maximum likelihood estimation. 
m1FULLML = lmer(difference_HR_log ~ duration*area +  wind_speed + temp + incubation_stage + (1|eider:bear_event) + (1|eider:bear_event:sample), data = NaturalPredation, REML=FALSE) 
summary(m1FULLML)
# interaction insignificant, remove from next model 
#Since slope is not significant, proceed with backward stepwise elimination removing non-significant terms, and using Maximum likelihood estimation. 
m1FULLML = lmer(difference_HR_log ~ duration +area +  wind_speed + temp + incubation_stage + (1|eider:bear_event) + (1|eider:bear_event:sample), data = NaturalPredation, REML=FALSE) 
summary(m1FULLML)

m1aML = lmer(difference_HR_log ~ bear_distance + duration_viewshed_scale + wind_speed + temp +  incubation_stage + (1|eider:bear_event) + (1|eider:bear_event:sample), data = NaturalPredation, REML=FALSE) 
summary(m1aML)
# bear distance most insignificant, remove from next model

m1bML = lmer(difference_HR_log ~  duration_viewshed_scale + wind_speed + temp +  incubation_stage + (1|eider:bear_event) + (1|eider:bear_event:sample), data = NaturalPredation, REML=FALSE) 
summary(m1bML)
#wind speed most insignificant, remove from next model 

m1cML = lmer(difference_HR_log ~ duration_viewshed_scale + temp +  incubation_stage + (1|eider:bear_event) + (1|eider:bear_event:sample), data = NaturalPredation, REML=FALSE) 
summary(m1cML)

m1cML = lmer(difference_HR_log ~ duration * area + temp +  incubation_stage + (1|eider:bear_event) + (1|eider:bear_event:sample), data = NaturalPredation, REML=FALSE) 
summary(m1cML)
#temp and incubation significant at 5% and duration significant at 5%
# was looking for a model where all parameters significant at the 0.1 level. 

#null model with only random effects
m1nullML = lmer(difference_HR_log ~ (1|eider:bear_event) + (1|eider:bear_event:sample), data = NaturalPredation, REML=FALSE) 
summary(m1nullML)

#AICc model selection with all models included in backward elimination process, and null model
model.sel(m1FULLML, m1aML, m1bML, m1cML, m1nullML)
# m1c and m1b equal 
m1dML = lmer(difference_HR_log ~ AVG.Baseline + duration_viewshed_scale + temp +  incubation_stage + (1|eider:bear_event) + (1|eider:bear_event:sample), data = NaturalPredation, REML=FALSE) 
summary(m1cML)
#To do: model selection of two models (AICc). magnitude and predHR <- (baseline + …) 
m1dML = lmer(predation_HR ~ AVG.Baseline + duration_viewshed_scale + temp +  incubation_stage + (1|eider:bear_event) + (1|eider:bear_event:sample), data = NaturalPredation, REML=FALSE) 
summary(m1cML)

#response variable differs between models
model.sel(m1cML, m1dML)

#refit m1c to REML
m1c = lmer(difference_HR_log ~ duration_viewshed_scale + temp +  incubation_stage + (1|eider:bear_event) + (1|eider:bear_event:sample), data = NaturalPredation, REML=TRUE) 
summary(m1c)
anova(m1c)
r.squaredGLMM(m1c)
qqnorm(residuals(m1c))
qqline(residuals(m1c))
hist(residuals(m1c))
boxplot(residuals(m1c))
shapiro.test(residuals(m1c))
vif(m1c)


# Plot predicted values 

mydfDur <- ggpredict(m1c,terms="duration_viewshed_scale")

ggplot(mydfDur, aes(x = x, y = predicted)) +
  geom_smooth(method="lm", 
              method.args = list(family = "gaussian"), colour = "black", se = FALSE, size=3) +
  geom_point(data=NaturalPredation,aes(x=duration_viewshed_scale, y=difference_HR_log), size=2, alpha=0.7, colour="black")+
  labs(y = "Log ∆HR (beats/10s)",
       x= bquote(bold("Eider Ratio of Exposure Risk to Polar Bear"~(min/m^2)))) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), linetype="dashed", alpha = 0.0, colour="black") +
  scale_x_continuous(breaks = seq(0, 220, by = 20), limits = c(-0, 220))+
  theme_classic() + theme(axis.text.x=element_text(size=12),axis.title=element_text(size=14))+
  theme(axis.text.x=element_text(size=26, family="Times New Roman"),
        axis.title=element_text(size=30,face = "bold", family="Times New Roman"))+ 
  theme(legend.title = element_text(size = 30, face="bold", family="Times New Roman"), 
        legend.text = element_text(size = 18))+
  theme(axis.text.y=element_text(size=26, family="Times New Roman")) + 
  theme(axis.line = element_line(size = 2, colour = "black", linetype=1))+
  scale_y_continuous(breaks = seq(0.4, 1.8, by = 0.2), limits = c(0.4, 1.8), expand=c(0,0))

#2) plot HR_difference to incubation stage 
mydfInc <- ggpredict(m1c,terms="incubation_stage")
plot(mydfInc)

ggplot(mydfInc, aes(x = x, y = predicted)) +
  geom_smooth(method="lm", 
              method.args = list(family = "gaussian"), colour = "black", se = FALSE, size=3, fullrange=TRUE) +
  geom_point(data=NaturalPredation,aes(x=incubation_stage, y=difference_HR_log), size=2, alpha=0.7, colour="black")+
  labs(y = "Log ∆HR (beats/10s)",
       x="Incubation Stage (days)") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), linetype="dashed", alpha = 0.0, colour="black") +
  theme_classic() + theme(axis.text.x=element_text(size=12),axis.title=element_text(size=14))+
  theme(axis.text.x=element_text(size=26, family="Times New Roman"),
        axis.title=element_text(size=30,face = "bold", family="Times New Roman"))+ 
  theme(legend.title = element_text(size = 30, face="bold", family="Times New Roman"), 
        legend.text = element_text(size = 18))+
  theme(axis.text.y=element_text(size=26,family="Times New Roman")) + 
  theme(axis.line = element_line(size = 2, colour = "black", linetype=1)) +
  scale_y_continuous(breaks = seq(0.4, 1.8, by = 0.2), limits = c(0.4, 1.8), expand=c(0,0))+
  scale_x_continuous(breaks = seq(4, 20, by = 1), limits = c(4, 19))

#3) plot HR_difference to temp 
mydfT <- ggpredict(m1c,terms="temp")
plot(mydfT)

ggplot(mydfT, aes(x = x, y = predicted)) +
  geom_smooth(method="lm", 
              method.args = list(family = "gaussian"), colour = "black", se = FALSE, size=3) +
  geom_point(data=NaturalPredation,aes(x=temp, y=difference_HR_log), size=2, alpha=0.7, colour="black")+
  labs(y = "Log ∆HR (beats/10s)",
       x="Air Temperature (degrees Celsius)") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), linetype="dashed", alpha = 0.0, colour="black") +
  scale_x_continuous(breaks = seq(0, 10, by = 1), limits = c(0, 10))+
  theme_classic() + theme(axis.text.x=element_text(size=12),axis.title=element_text(size=14))+
  theme(axis.text.x=element_text(size=30,family="Times New Roman"),
        axis.title=element_text(size=30,face = "bold", family="Times New Roman"))+ 
  theme(legend.title = element_text(size = 26, face="bold",family="Times New Roman"), 
        legend.text = element_text(size = 26, family="Times New Roman"))+
  theme(axis.text.y=element_text(size=30,family="Times New Roman")) + 
  theme(axis.line = element_line(size = 2, colour = "black", linetype=1)) +
  scale_y_continuous(breaks = seq(0.4, 1.8, by = 0.2), limits = c(0.4, 1.8), expand=c(0,0))


#- Compare baseline heart rate to predator-induced heart rate (LMM)

HRType <- read.csv(file.choose())

HRType <-HRType %>%
  rename(eider = Recorder.ID,
         interval = Interval.Number,
         bear_event = Bear.Event.b,
         heart_rate = Total.Count,
         sample = Sample.Number,
         HR_type = Heart.Rate.Type)

HRType$eider<-as.factor(HRType$eider)  
HRType$interval<-as.factor(HRType$interval)
HRType$bear_event<-as.factor(HRType$bear_event)
HRType$sample<-as.factor(HRType$sample)

HR <- lmer(heart_rate ~ HR_type + (1|eider:bear_event) + (1|eider:bear_event:sample), data = HRType, REML=TRUE) 
summary(HR)
r.squaredGLMM(HR)
#heart rate to predator is  significant higher than to baseline
#posthoc testing to confirm direction 
lsmeans(HR, pairwise ~ HR_type, adjust= "tukey", level=0.95)

#Summary statistics
RTypeB <- HRType %>%
  filter(HR_type == "Baseline") 
summary(HRTypeB$heart_rate)


HRTypeP <- HRType %>%
  filter(HR_type == "Predator") 
summary(HRTypeP$heart_rate)


#Figure comparing baseline to bear-induced heart rate

mydfHR <- ggpredict(HR, terms = "HR_type")

ggplot(HRType, aes(x = HR_type, y = heart_rate)) +
  geom_boxplot() + 
  theme_classic() +
  theme(axis.text.x=element_text(size=12),axis.title=element_text(size=14))+
  theme(axis.text.x=element_text(size=18),
        axis.title=element_text(size=24,face = "bold"))+ 
  theme(legend.title = element_text(size = 24, face="bold"), 
        legend.text = element_text(size = 18))+
  theme(axis.text.y=element_text(size=18))+
  labs(y = "Heart Rate",
       x="Response Type") + 
  theme(text = element_text(family = "Times New Roman")) + geom_point(data=mydfHR, aes(x = x, y =  predicted), position= position_nudge(x=0), colour='red')+
  geom_segment(data=mydfHR, aes(x=x, xend=x, y=conf.low, yend=conf.high), colour='red') +
scale_x_discrete(labels = c('Baseline','Polar Bear'))

#- Eider heart rate over time (intervalB) within a bear event
#Change sample interval to an ordinal variable 
NaturalPredation$intervalB <- ordered(NaturalPredation$intervalB, levels = c("1", "2", "3","4", "5", "6", "7", "8", "9", "10","11", "12", "13", "14", "15","16", "17", "18", "19", "20", "21", "22","23", "24", "25", "26", "27","28", "29", "30", "31", "32", "33", "34", "35","36"))

modelC = lmer(predation_HR ~ intervalB + (1|eider:bear_event), 
              data = NaturalPredation)
summary(modelC)
anova(modelC)

#create a figure 

mydfintervalB <- ggpredict(modelC, terms = "intervalB")

ggplot(NaturalPredation, aes(x = intervalB, y = predation_HR)) +
  geom_boxplot() + 
  geom_point(data=mydfintervalB, aes(x = x, y =  predicted), position= position_nudge(x=0), colour='red')+
  geom_segment(data=mydfintervalB, aes(x=x, xend=x, y=conf.low, yend=conf.high), colour='red')+ 
  theme_classic() +
  theme(axis.text.x=element_text(size=12),axis.title=element_text(size=14))+
  theme(axis.text.x=element_text(size=18),
        axis.title=element_text(size=24,face = "bold"))+ 
  theme(legend.title = element_text(size = 24, face="bold"), 
        legend.text = element_text(size = 18))+
  theme(axis.text.y=element_text(size=18))+
  labs(y = "Heart Rate (beats/10s)",
       x="Sample Interval") +
  theme(text = element_text(family = "Times New Roman"))+
  theme(strip.text = element_text(
    size = 18, face="bold"))

#Check for correlation between area of visibility and exposure duration

Exposure <- read.csv(file.choose())

Exposure <-Exposure %>%
  rename(eider = Recorder.ID,
         bear_event = Bear.Event.b,
         duration = Duration.Bear.Viewshed..min.,
         area = Viewshed.Area..m2.)

shapiro.test(Exposure$duration) # => p < 0.05 data are not normal
shapiro.test(Exposure$area) # => p < 0.05 data are not normal

cor.test(Exposure$duration, Exposure$area, method = "spearman") #duration and area are correlated 
#Spearman is more appropriate for measurements taken from ordinal scales.

ggplot(NaturalPredation, aes(x = area, y = duration, color=duration_viewshed_area)) +
  geom_point(size=3) +scale_color_gradientn(colours = rainbow(5)) + 
  theme_classic() + 
  labs(y = "Exposure Duration (min)",
       x=bquote(bold("Viewshed Area"~(m^2))), color="Ratio of risk") + theme(axis.text.x=element_text(size=12),axis.title=element_text(size=14))+
  theme(axis.text.x=element_text(size=26, family="Times New Roman"),
        axis.title=element_text(size=30,face = "bold", family="Times New Roman"))+ 
  theme(legend.title = element_text(size = 30, face="bold", family="Times New Roman"), 
        legend.text = element_text(size = 18))+
  theme(axis.text.y=element_text(size=26, family="Times New Roman")) + 
  theme(axis.line = element_line(size = 2, colour = "black", linetype=1))

#####Supplemental figures 

# Plot predicted values 

mydfDur <- ggpredict(m1c,terms="duration_viewshed_scale")
mydfDur$predicted_back <- (10^(mydfDur$predicted)-10)
mydfDur$conf.low_back <- (10^(mydfDur$conf.low)-10)
mydfDur$conf.high_back <- (10^(mydfDur$conf.high)-10)

ggplot()  +
  labs(y = "Predicted ∆HR (beats/10s)",
       x= bquote(bold("Eider Ratio of Exposure Risk to Polar Bear"~(min/m^2)))) +
  scale_x_continuous(breaks = seq(0, 220, by = 20), limits = c(-0, 220))+
  theme_classic() + theme(axis.text.x=element_text(size=12),axis.title=element_text(size=14))+
  theme(axis.text.x=element_text(size=26, family="Times New Roman"),
        axis.title=element_text(size=30,face = "bold", family="Times New Roman"))+ 
  theme(legend.title = element_text(size = 30, face="bold", family="Times New Roman"), 
        legend.text = element_text(size = 18))+
  theme(axis.text.y=element_text(size=26, family="Times New Roman")) + 
  theme(axis.line = element_line(size = 2, colour = "black", linetype=1)) +
  geom_point(data=mydfDur, aes(x = x, y =  predicted_back), position= position_nudge(x=0), colour='black', size=4)+
  scale_y_continuous(breaks = seq(-4, 4, by = 2), limits = c(-4, 5), expand=c(0,0))


#2) plot HR_difference to incubation stage 
mydfInc <- ggpredict(m1c,terms="incubation_stage")
mydfInc$predicted_back <- (10^(mydfInc$predicted)-10)
mydfInc$conf.low_back <- (10^(mydfInc$conf.low)-10)
mydfInc$conf.high_back <- (10^(mydfInc$conf.high)-10)
plot(mydfInc)

ggplot() +
  labs(y = "Predicted ∆HR (beats/10s)",
       x="Incubation Stage (days)") +
  theme_classic() + theme(axis.text.x=element_text(size=12),axis.title=element_text(size=14))+
  theme(axis.text.x=element_text(size=26, family="Times New Roman"),
        axis.title=element_text(size=30,face = "bold", family="Times New Roman"))+ 
  theme(legend.title = element_text(size = 30, face="bold", family="Times New Roman"), 
        legend.text = element_text(size = 18))+
  theme(axis.text.y=element_text(size=26,family="Times New Roman")) + 
  theme(axis.line = element_line(size = 2, colour = "black", linetype=1)) +
  scale_x_continuous(breaks = seq(4, 20, by = 1), limits = c(4, 19))+
  geom_point(data=mydfInc, aes(x = x, y =  predicted_back), position= position_nudge(x=0), colour='black', size=4) +
scale_y_continuous(breaks = seq(-2, 10, by = 2), limits = c(-3, 11), expand=c(0,0))

#3) plot HR_difference to temp 
mydfT <- ggpredict(m1c,terms="temp")
mydfT$predicted_back <- (10^(mydfT$predicted)-10)
mydfT$conf.low_back <- (10^(mydfT$conf.low)-10)
mydfT$conf.high_back <- (10^(mydfT$conf.high)-10)
plot(mydfT)

ggplot() +
  labs(y = "Predicted ∆HR (beats/10s)",
       x="Air Temperature (degrees Celsius)") +
  scale_x_continuous(breaks = seq(0, 10, by = 1), limits = c(0, 10))+
  theme_classic() + theme(axis.text.x=element_text(size=12),axis.title=element_text(size=14))+
  theme(axis.text.x=element_text(size=30,family="Times New Roman"),
        axis.title=element_text(size=30,face = "bold", family="Times New Roman"))+ 
  theme(legend.title = element_text(size = 26, face="bold",family="Times New Roman"), 
        legend.text = element_text(size = 26, family="Times New Roman"))+
  theme(axis.text.y=element_text(size=30,family="Times New Roman")) + 
  theme(axis.line = element_line(size = 2, colour = "black", linetype=1)) +
  geom_point(data=mydfT, aes(x = x, y =  predicted_back), position= position_nudge(x=0), colour='black', size=4)+
  scale_y_continuous(breaks = seq(-2, 6, by = 2), limits = c(-3, 7), expand=c(0,0))