## LOADING PACKAGES
library(ggplot2)
#install.packages("glmmTMB")
library(glmmTMB)
#install.packages("pscl")
library(pscl)
#install.packages("bbmle")
library(bbmle)
#install.packages("dunn.test")
library(dunn.test)
#install.packages("lubridate")
library(lubridate)
#install.packages("activity")
library(activity)
#install.packages("ggeffects")
library(ggeffects)
#install.packages("sjPlot")
library(sjPlot)
#install.packages("DHARMa")
library(DHARMa)
#install.packages("SciViews")
library(SciViews)

# LOADING DATA
data <- read.table("/home/peterxps/Dropbox/papers/Clamsen_Kori_Bustard_movement_patterns/Kori_GPS_CJedited_PSR_verified4.csv", sep=";", header=TRUE)

str(data)
head(data)

# convert month to factor
data$Month <- as.factor(as.integer(data$Month))

# sample sizes pr month, sex and day-category
table(data$Month)
table(data$CollarID, data$Sex)
table(data$Mean_cat)

# round distance pr hr and make histogram
data$Dist_pr_hr <- round(data$Dist_pr_hr, 0)
hist(data$Dist_pr_hr, breaks=50)

# log distances
data$logdist <- log(data$Dist_pr_hr+1)
hist(data$logdist)

# Visualise each collar
levels(data$CollarID)
hist(data$logdist[data$CollarID == "1279"])
hist(data$logdist[data$CollarID == "1280"])
hist(data$logdist[data$CollarID == "1281"])
hist(data$logdist[data$CollarID == "1282"])
hist(data$logdist[data$CollarID == "631"])
hist(data$logdist[data$CollarID == "632"])
hist(data$logdist[data$CollarID == "633"])
hist(data$logdist[data$CollarID == "634"])

# Table S1
tapply(data$Dist_pr_hr, data$CollarID, mean, na.rm=T)
tapply(data$Dist_pr_hr, data$CollarID, sd, na.rm=T)
tapply(data$Dist_pr_hr, data$CollarID, max, na.rm=T)
tapply(data$Dist_pr_hr, data$CollarID, min, na.rm=T)
tapply(data$Dist_pr_hr, data$Mean_cat, mean, na.rm=T)
tapply(data$Dist_pr_hr, data$Mean_cat, sd, na.rm=T)
tapply(data$Dist_pr_hr, data$Month, mean, na.rm=T)

# Overall mean and SD
mean(data$Dist_pr_hr, na.rm=T)
sqrt(var(data$Dist_pr_hr, na.rm=T))


# hrs since last reading
table(data$Hrs_since_last)
median(data$Hrs_since_last, na.rm = T)
hist(data$Hrs_since_last[data$Hrs_since_last<13])

# set morning as reference level
data$Mean_cat <- relevel(data$Mean_cat, ref = "1_Morning")

# make subset data2, where time between readings is 8 hrs or less
data2 <- subset(data, Hrs_since_last < 9)
median(data2$Hrs_since_last)



# Overall mean and SD (data2)
mean(data2$Dist_pr_hr, na.rm=T)
sqrt(var(data2$Dist_pr_hr, na.rm=T))

# Table S1 (data2)
tapply(data2$Dist_pr_hr, data2$Mean_cat, mean, na.rm=T)
tapply(data2$Dist_pr_hr, data2$Mean_cat, sd, na.rm=T)
tapply(data2$Dist_pr_hr, data2$Month, mean, na.rm=T)







#Table 1
data2$Mean_cat <- relevel(data2$Mean_cat, ref = "1_Morning")

mod2 <- glmmTMB(Dist_pr_hr ~ Mean_cat + (1|CollarID), data = data2, family = nbinom1)
summary(mod2)
mod_simres <- simulateResiduals(mod2, quantreg=T)
plot(mod_simres)




#Table S2 - includinig sex
mod2s2 <- glmmTMB(Dist_pr_hr ~ Mean_cat + Sex + (1|CollarID), data = data2, family = nbinom1)
summary(mod2s2)
mod_simress2 <- simulateResiduals(mod2s2, quantreg=T)
plot(mod_simress2)

#Table S3 - including month
mod2s3 <- glmmTMB(Dist_pr_hr ~ Mean_cat + Month + (1|CollarID), data = data2, family = nbinom1)
summary(mod2s3)
mod_simress3 <- simulateResiduals(mod2s3, quantreg=T)
plot(mod_simress3)


# NOTE slight deviations from model assumptions

# Three figures respectively - no effect of including sex or month
p1 <- plot_model(mod2, 
                 type = "pred", 
                 terms = c("Mean_cat"), 
                 dot.size = 3, 
                 line.size = 2.0, 
                 width = 0.5, 
                 axis.title = c("Period of the day", "Distance moved per hour (m)"),
                 title = "",
                 show.data = FALSE,
                 show.legend = TRUE
) + theme_sjplot()
p1

p2 <- plot_model(mod2s2, 
                 type = "pred", 
                 terms = c("Mean_cat"), 
                 dot.size = 3, 
                 line.size = 2.0, 
                 width = 0.5, 
                 axis.title = c("Period of the day", "Distance moved per hour (m)"),
                 title = "",
                 show.data = FALSE,
                 show.legend = TRUE
) + theme_sjplot()
p2

p3 <- plot_model(mod2s3, 
                 type = "pred", 
                 terms = c("Mean_cat"), 
                 dot.size = 3, 
                 line.size = 2.0, 
                 width = 0.5, 
                 axis.title = c("Period of the day", "Distance moved per hour (m)"),
                 title = "",
                 show.data = FALSE,
                 show.legend = TRUE
) + theme_sjplot()
p3