####################################################
# ERT Movement Analysis
#
# Author: Dr Vicky Boult
# Date: March 2023
####################################################

# Load relevant packages
library(dplyr)
library(ggplot2)

# Read in data
kkr <- read.csv("C:/Users/ku912605/Dropbox/Elephant_Reintergration_Trust/ERT Movement/KKR_movement_metrics.csv", header = T)
spgr <- read.csv("C:/Users/ku912605/Dropbox/Elephant_Reintergration_Trust/ERT Movement/SPGR_movement_metrics.csv", header = T)
hr <- read.csv("C:/Users/ku912605/Dropbox/Elephant_Reintergration_Trust/ERT Movement/KKR_SPGR_homerange.csv", header = T)

kkr <- read.csv("E:/Downloads/OD Personal/Desktop/KKR_movement_metrics.csv", header = T)
spgr <- read.csv("E:/Downloads/OD Personal/Desktop/SPGR_movement_metrics.csv", header = T)

#----------------------------------------------------
# Pre-processing
#----------------------------------------------------
# Combine KKR and SPGR data frames for easy processing
all_move <- rbind.data.frame(kkr, spgr)

# Date & time formatting 
all_move$Date <- as.Date(all_move$Date, format = "%d/%m/%Y")
all_move$datetime <- as.POSIXct(all_move$timestamp, format = "%d/%m/%Y %H:%M")

# Drop four dodgey rows (datetime format not working..?)
all_move <- all_move %>%
  filter(!row_number() %in% c(1883, 12176, 28895, 30324))

# Format phase as categorical
all_move$Phase <- as.factor(all_move$Phase)

# Calculate the daily total distance traveled
daily_dist <- all_move %>%
  group_by(animal, Date) %>%
  mutate(total_dist = sum(distance)) %>%
  filter(row_number()==1) %>%
  filter(!is.na(total_dist))

# Separate into KKR and SPGR for analysis
kkr_dd <- daily_dist %>%
  filter(animal == "kkr_herd_bull")
spgr_dd <- daily_dist %>%
  filter(animal == "spgr_herd_bull")

# Calculate distance traveled for each fix accounting for time between subsequent locations
time_dist <- all_move %>%
  mutate(diff = datetime - lag(datetime),
         diff_secs = as.numeric(diff, units = "secs"),
         dist_per_sec = distance / diff_secs)

# Assign times to "period" of day (i.e. day, night, dawn, dusk)
time_dist$time_new <- format(as.POSIXct(time_dist$datetime), format = "%H:%M")
# (These times can be adjusted as desired)
dawn_start <- format(as.POSIXct("05:00", format = "%H:%M"), format = "%H:%M")
dawn_end <- format(as.POSIXct("07:00", format = "%H:%M"), format = "%H:%M")
day_start <- format(as.POSIXct("12:00", format = "%H:%M"), format = "%H:%M")
day_end <- format(as.POSIXct("14:00", format = "%H:%M"), format = "%H:%M")
dusk_start <- format(as.POSIXct("18:00", format = "%H:%M"), format = "%H:%M")
dusk_end <- format(as.POSIXct("20:00", format = "%H:%M"), format = "%H:%M")
night_start <- format(as.POSIXct("22:00", format = "%H:%M"), format = "%H:%M")
night_end <- format(as.POSIXct("02:00", format = "%H:%M"), format = "%H:%M")
# Create new variable
time_dist$period <- NA
# Fill new variable with period of the day
for (i in c(1:nrow(time_dist))){
  if (time_dist$time_new[i] >= dawn_start && time_dist$time_new[i] < dawn_end) {time_dist$period[i] <- "dawn"} 
  else if (time_dist$time_new[i] >= day_start && time_dist$time_new[i] < day_end) {time_dist$period[i] <- "day"}
  else if (time_dist$time_new[i] >= dusk_start && time_dist$time_new[i] < dusk_end) {time_dist$period[i] <- "dusk"}
  else if (time_dist$time_new[i] >= night_start || time_dist$time_new[i] < night_end) {time_dist$period[i] <- "night"}}

# Separate into KKR and SPGR for analysis
kkr_td <- time_dist %>%
  filter(animal == "kkr_herd_bull")
spgr_td <- time_dist %>%
  filter(animal == "spgr_herd_bull")

# Separate home range by reserve
kkr_hr <- hr %>%
  filter(Reserve == "KKR")
spgr_hr <- hr %>%
  filter(Reserve == "SPGR")

#----------------------------------------------------
# 1. Phase comparison
#----------------------------------------------------
# SEASONALITY
# Two-way ANOVA to compare differences between phases, seasons and interaction

#------
# KKR
#------

# First, test for normality
shapiro.test(kkr_dd$total_dist) # Not normal
hist(kkr_dd$total_dist) # Left skewed - try log transform
shapiro.test(log(kkr_dd$total_dist)) # Success!

# Perform two-way ANOVA
aov_dd <- aov(log(total_dist) ~ Season * Phase, data = kkr_dd)
#plot(aov_dd) # Check for assumptions of normality and homoscedascity
summary(aov_dd)
TukeyHSD(aov_dd)

# Plot seasonal means per Phase
kkr_dd %>%
  mutate(Season = factor(Season, levels = c("Wet", "Fall", "Dry", "Spring"))) %>%
  ggplot() +
  geom_boxplot(aes(x = Season, y = log(total_dist), fill = Phase), lwd = 0.75) +
  scale_fill_manual(values = c("indianred2", "dodgerblue")) +
  xlab("Season") + ylab("Log total distance (m)") +
  theme_bw() +
  theme(axis.title = element_text(size = 12, face = "bold"),
        axis.text = element_text(size = 11),
        legend.title = element_text(size = 12, face = "bold"),
        legend.text = element_text(size = 11))

#------
# SPGR
#------

# First, test for normality
shapiro.test(spgr_dd$total_dist) # Not normal
hist(spgr_dd$total_dist) # Left skewed - try log transform
shapiro.test(log(spgr_dd$total_dist)) # Close enough (qqplot of residuals looks good)

# Perform two-way ANOVA
aov_dd <- aov(log(total_dist) ~ Season * Phase, data = spgr_dd)
#plot(aov_dd) # Check for assumptions of normality and homoscedascity
summary(aov_dd)
TukeyHSD(aov_dd)

# Plot seasonal means per Phase
spgr_dd %>%
  mutate(Season = factor(Season, levels = c("Wet", "Fall", "Dry", "Spring"))) %>%
  ggplot() +
  geom_boxplot(aes(x = Season, y = log(total_dist), fill = Phase), lwd = 0.75) +
  scale_fill_manual(values = c("indianred2", "dodgerblue")) +
  xlab("Season") + ylab("Log total distance (m)") +
  theme_bw() +
  theme(axis.title = element_text(size = 12, face = "bold"),
        axis.text = element_text(size = 11),
        legend.title = element_text(size = 12, face = "bold"),
        legend.text = element_text(size = 11))

#----------------------
# TIME OF DAY
#----------------------
# Not normally distributed so can't do a two-way ANOVA

#------
# KKR
#------

# Calculate mean distance per second for each period
kkr_td_mean <- kkr_td %>%
  filter(!is.na(dist_per_sec), !is.na(period)) %>%
  group_by(Phase, period) %>%
  summarise(mean_dist = mean(dist_per_sec, na.rm = T))

# Build dataframe for chi-squared
kkr_td_chi <- cbind.data.frame(kkr_td_mean[1:4,2], kkr_td_mean[1:4,3], kkr_td_mean[5:8,3])
colnames(kkr_td_chi) <- c("period", "phase1", "phase2")

# Chi-squared
chisq.test(kkr_td_chi[,2:3]) # Not significantly different

# Plot distance travelled per Phase
kkr_td %>%
  filter(!is.na(period)) %>%
  ggplot() +
  geom_boxplot(aes(x = period, y = log(dist_per_sec), fill = Phase), lwd = 0.75) +
  scale_fill_manual(values = c("indianred2", "dodgerblue")) +
  scale_x_discrete(labels = c("Dawn", "Day", "Dusk", "Night")) +
  xlab("Season") + ylab("Log distance travelled (m/s)") +
  theme_bw() +
  theme(axis.title = element_text(size = 12, face = "bold"),
        axis.text = element_text(size = 11),
        legend.title = element_text(size = 12, face = "bold"),
        legend.text = element_text(size = 11))

#------
# SPGR
#------

# Calculate mean distance per second for each period
spgr_td_mean <- spgr_td %>%
  filter(!is.na(dist_per_sec), !is.na(period)) %>%
  group_by(Phase, period) %>%
  summarise(mean_dist = mean(dist_per_sec, na.rm = T))

# Build dataframe for chi-squared
spgr_td_chi <- cbind.data.frame(spgr_td_mean[1:4,2], spgr_td_mean[1:4,3], spgr_td_mean[5:8,3])
colnames(spgr_td_chi) <- c("period", "phase1", "phase2")

# Chi-squared
chisq.test(spgr_td_chi[,2:3]) # Not significantly different

# Plot distance travelled per Phase
spgr_td %>%
  filter(!is.na(period)) %>%
  ggplot() +
  geom_boxplot(aes(x = period, y = log(dist_per_sec), fill = Phase), lwd = 0.75) +
  scale_fill_manual(values = c("indianred2", "dodgerblue")) +
  scale_x_discrete(labels = c("Dawn", "Day", "Dusk", "Night")) +
  xlab("Season") + ylab("Log distance travelled (m/s)") +
  theme_bw() +
  theme(axis.title = element_text(size = 12, face = "bold"),
        axis.text = element_text(size = 11),
        legend.title = element_text(size = 12, face = "bold"),
        legend.text = element_text(size = 11))

#----------------------
# HOME RANGE
#----------------------

#------
# KKR
#------

# Calculate mean home range for each Phase * Season combination
kkr_hr <- kkr_hr %>%
  group_by(Phase, Season) %>%
  summarise(mean_hr = mean(Area))
kkr_hr$Phase <- as.character(kkr_hr$Phase)

# Build chi-squared data frame
kkr_hr_chi <- cbind.data.frame(kkr_hr[1:4,2], kkr_hr[1:4,3], kkr_hr[5:8,3])
colnames(kkr_hr_chi) <- c("season", "phase1", "phase2")

# Chi-squared 
chisq.test(kkr_hr_chi[,2:3])

# Plot seasonal home range size per phase
kkr_hr %>%
  mutate(Season = factor(Season, levels = c("Wet", "Fall", "Dry", "Spring"))) %>%
  ggplot(aes(x = Season, y = mean_hr, group = Phase, fill = Phase)) +
  geom_bar(stat = "summary", fun = "mean", width = 0.5, position = position_dodge()) +
  scale_fill_manual(values = c("indianred2", "dodgerblue")) +
  xlab("Season") + ylab("Home range size") +
  theme_bw() +
  theme(axis.title = element_text(size = 12, face = "bold"),
        axis.text = element_text(size = 11),
        legend.title = element_text(size = 12, face = "bold"),
        legend.text = element_text(size = 11))

#------
# SPGR
#------

# Calculate mean home range for each Phase * Season combination
spgr_hr <- spgr_hr %>%
  group_by(Phase, Season) %>%
  summarise(mean_hr = mean(Area))
spgr_hr$Phase <- as.character(spgr_hr$Phase)

# Build chi-squared data frame
spgr_hr_chi <- cbind.data.frame(spgr_hr[1:4,2], spgr_hr[1:4,3], spgr_hr[5:8,3])
colnames(spgr_hr_chi) <- c("season", "phase1", "phase2")

# Chi-squared 
chisq.test(spgr_hr_chi[,2:3])

# Plot seasonal home range size per phase
spgr_hr %>%
  mutate(Season = factor(Season, levels = c("Wet", "Fall", "Dry", "Spring"))) %>%
  ggplot(aes(x = Season, y = mean_hr, group = Phase, fill = Phase)) +
  geom_bar(stat = "summary", fun = "mean", width = 0.5, position = position_dodge()) +
  scale_fill_manual(values = c("indianred2", "dodgerblue")) +
  xlab("Season") + ylab("Home range size") +
  theme_bw() +
  theme(axis.title = element_text(size = 12, face = "bold"),
        axis.text = element_text(size = 11),
        legend.title = element_text(size = 12, face = "bold"),
        legend.text = element_text(size = 11))

#----------------------------------------------------
# 2. Match wild eles?
#----------------------------------------------------

#----------------------
# SEASONALITY
#----------------------

#------
# KKR
#------

# Split data frame into Phase 1 and 2
kkr_dd_1 <- kkr_dd %>%
  filter(Phase == 1)
kkr_dd_2 <- kkr_dd %>%
  filter(Phase == 2)

# PHASE 1

# Test for normality
shapiro.test(kkr_dd_1$total_dist) # Not normal
hist(kkr_dd_1$total_dist) # Left skewed - try log transform
shapiro.test(log(kkr_dd_1$total_dist)) # Success!

# ANOVA log daily distance by season
season_kkr_1 <- aov(log(total_dist) ~ Season, data = kkr_dd_1) # Log total daily distance to normalise
#plot(season_kkr_1)
summary(season_kkr_1) # Significant
TukeyHSD(season_kkr_1) # Post-hoc to identify where differences lie

# Plot daily distance per season - ensure it aligns with Tukey
kkr_dd_1 %>%
  mutate(Season = factor(Season, levels = c("Wet", "Fall", "Dry", "Spring"))) %>%
  ggplot() +
  geom_boxplot(aes(x = Season, y = log(total_dist)), fill = "indianred2", lwd = 0.75) +
  #scale_fill_manual(values = c("indianred2", "dodgerblue")) +
  xlab("Season") + ylab("Log total distance (m)") +
  theme_bw() +
  theme(axis.title = element_text(size = 12, face = "bold"),
        axis.text = element_text(size = 11),
        legend.title = element_text(size = 12, face = "bold"),
        legend.text = element_text(size = 11))

# PHASE 2

# Test for normality
shapiro.test(kkr_dd_2$total_dist) # Not normal
hist(kkr_dd_2$total_dist) # Left skewed - try log transform
shapiro.test(log(kkr_dd_2$total_dist)) # Success!

# ANOVA log daily distance by season
season_kkr_2 <- aov(log(total_dist) ~ Season, data = kkr_dd_2) # Log total daily distance to normalise
#plot(season_kkr_2)
summary(season_kkr_2) # Significant
TukeyHSD(season_kkr_2) # Post-hoc to identify where differences lie

# Plot daily distance per season - ensure it aligns with Tukey
kkr_dd_2 %>%
  mutate(Season = factor(Season, levels = c("Wet", "Fall", "Dry", "Spring"))) %>%
  ggplot() +
  geom_boxplot(aes(x = Season, y = log(total_dist)), fill = "indianred2", lwd = 0.75) +
  #scale_fill_manual(values = c("indianred2", "dodgerblue")) +
  xlab("Season") + ylab("Log total distance (m)") +
  theme_bw() +
  theme(axis.title = element_text(size = 12, face = "bold"),
        axis.text = element_text(size = 11),
        legend.title = element_text(size = 12, face = "bold"),
        legend.text = element_text(size = 11))

#------
# SPGR
#------

# Split data frame into Phase 1 and 2
spgr_dd_1 <- spgr_dd %>%
  filter(Phase == 1)
spgr_dd_2 <- spgr_dd %>%
  filter(Phase == 2)

# PHASE 1

# Test for normality
shapiro.test(spgr_dd_1$total_dist) # Not normal
hist(spgr_dd_1$total_dist) # Left skewed - try log transform
shapiro.test(log(spgr_dd_1$total_dist)) # Close enough!

# ANOVA log daily distance by season
season_spgr_1 <- aov(log(total_dist) ~ Season, data = spgr_dd_1) # Log total daily distance to normalise
#plot(season_spgr_1)
summary(season_spgr_1) # Significant
TukeyHSD(season_spgr_1) # Post-hoc to identify where differences lie

# Plot daily distance per season - ensure it aligns with Tukey
spgr_dd_1 %>%
  mutate(Season = factor(Season, levels = c("Wet", "Fall", "Dry", "Spring"))) %>%
  ggplot() +
  geom_boxplot(aes(x = Season, y = log(total_dist)), fill = "indianred2", lwd = 0.75) +
  #scale_fill_manual(values = c("indianred2", "dodgerblue")) +
  xlab("Season") + ylab("Log total distance (m)") +
  theme_bw() +
  theme(axis.title = element_text(size = 12, face = "bold"),
        axis.text = element_text(size = 11),
        legend.title = element_text(size = 12, face = "bold"),
        legend.text = element_text(size = 11))

# PHASE 2

# Test for normality
shapiro.test(spgr_dd_2$total_dist) # Not normal
hist(spgr_dd_2$total_dist) # Left skewed - try log transform
shapiro.test(log(spgr_dd_2$total_dist)) # Success!

# ANOVA log daily distance by season
season_spgr_2 <- aov(log(total_dist) ~ Season, data = spgr_dd_2) # Log total daily distance to normalise
#plot(season_spgr_2)
summary(season_spgr_2) # Significant
TukeyHSD(season_spgr_2) # Post-hoc to identify where differences lie

# Plot daily distance per season - ensure it aligns with Tukey
spgr_dd_2 %>%
  mutate(Season = factor(Season, levels = c("Wet", "Fall", "Dry", "Spring"))) %>%
  ggplot() +
  geom_boxplot(aes(x = Season, y = log(total_dist)), fill = "indianred2", lwd = 0.75) +
  #scale_fill_manual(values = c("indianred2", "dodgerblue")) +
  xlab("Season") + ylab("Log total distance (m)") +
  theme_bw() +
  theme(axis.title = element_text(size = 12, face = "bold"),
        axis.text = element_text(size = 11),
        legend.title = element_text(size = 12, face = "bold"),
        legend.text = element_text(size = 11))

#----------------------
# TIME OF DAY
#----------------------

#------
# KKR
#------

# Split data frame into Phases 1 and 2
kkr_td_1 <- kkr_td %>%
  filter(Phase == 1, !is.na(period))
kkr_td_2 <- kkr_td %>%
  filter(Phase == 2, !is.na(period))

# PHASE 1

# Check for normality
shapiro.test(kkr_td_1$dist_per_sec) # Not normal
hist(kkr_td_1$dist_per_sec) # Heavily right skewed
shapiro.test(log(kkr_td_1$dist_per_sec)) 
# Tried a bunch of transformations, but no luck, use non-parametric alternative

install.packages("FSA")
library(FSA)

# Kruskal Wallis
kruskal.test(dist_per_sec ~ period, data = kkr_td_1) # Significant!
# Use Dunn post hoc to identify where significance lies
# Bonferroni's correction to account for multiple pairwise comparison 
dunnTest(dist_per_sec ~ period, data = kkr_td_1, method = "bonferroni")

# Plot distance per period - ensure it aligns with post-hoc
kkr_td_1 %>%
  filter(!is.na(period)) %>%
  ggplot() +
  geom_boxplot(aes(x = period, y = log(dist_per_sec)), fill = "indianred2", lwd = 0.75) +
  scale_x_discrete(labels = c("Dawn", "Day", "Dusk", "Night")) +
  xlab("Season") + ylab("Log distance travelled (m/s)") +
  theme_bw() +
  theme(axis.title = element_text(size = 12, face = "bold"),
        axis.text = element_text(size = 11),
        legend.title = element_text(size = 12, face = "bold"),
        legend.text = element_text(size = 11))

# PHASE 2 

# Check for normality
shapiro.test(kkr_td_2$dist_per_sec) # Not normal
hist(kkr_td_2$dist_per_sec) # Heavily right skewed
shapiro.test(log(kkr_td_2$dist_per_sec)) 
# Tried a bunch of transformations, but no luck, use non-parametric alternative


dunnTest.(FSA)

# Kruskal Wallis
kruskal.test(dist_per_sec ~ period, data = kkr_td_2) # Significant!
# Use Dunn post hoc to identify where significance lies
# Bonferroni's correction to account for multiple pairwise comparison 
dunnTest(dist_per_sec ~ period, data = kkr_td_2, method = "bonferroni")

# Plot distance per period - ensure it aligns with post-hoc
kkr_td_2 %>%
  filter(!is.na(period)) %>%
  ggplot() +
  geom_boxplot(aes(x = period, y = log(dist_per_sec)), fill = "indianred2", lwd = 0.75) +
  scale_x_discrete(labels = c("Dawn", "Day", "Dusk", "Night")) +
  xlab("Season") + ylab("Log distance travelled (m/s)") +
  theme_bw() +
  theme(axis.title = element_text(size = 12, face = "bold"),
        axis.text = element_text(size = 11),
        legend.title = element_text(size = 12, face = "bold"),
        legend.text = element_text(size = 11))

#------
# SPGR
#------

# Split data frame into Phases 1 and 2
spgr_td_1 <- spgr_td %>%
  filter(Phase == 1, !is.na(period))
spgr_td_2 <- spgr_td %>%
  filter(Phase == 2, !is.na(period))

# PHASE 1

# Check for normality
shapiro.test(spgr_td_1$dist_per_sec) # Not normal
hist(spgr_td_1$dist_per_sec) # Heavily right skewed
shapiro.test(log(spgr_td_1$dist_per_sec)) 
# Tried a bunch of transformations, but no luck, use non-parametric alternative

# Kruskal Wallis
kruskal.test(dist_per_sec ~ period, data = spgr_td_1) # Significant!
# Use Dunn post hoc to identify where significance lies
# Bonferroni's correction to account for multiple pairwise comparison 
dunnTest(dist_per_sec ~ period, data = spgr_td_1, method = "bonferroni")

# Plot distance per period - ensure it aligns with post-hoc
spgr_td_1 %>%
  filter(!is.na(period)) %>%
  ggplot() +
  geom_boxplot(aes(x = period, y = log(dist_per_sec)), fill = "indianred2", lwd = 0.75) +
  scale_x_discrete(labels = c("Dawn", "Day", "Dusk", "Night")) +
  xlab("Season") + ylab("Log distance travelled (m/s)") +
  theme_bw() +
  theme(axis.title = element_text(size = 12, face = "bold"),
        axis.text = element_text(size = 11),
        legend.title = element_text(size = 12, face = "bold"),
        legend.text = element_text(size = 11))

# PHASE 2 

# Check for normality
shapiro.test(spgr_td_2$dist_per_sec) # Not normal
hist(spgr_td_2$dist_per_sec) # Heavily right skewed
shapiro.test(log(spgr_td_2$dist_per_sec)) 
# Tried a bunch of transformations, but no luck, use non-parametric alternative

# Kruskal Wallis
kruskal.test(dist_per_sec ~ period, data = spgr_td_2) # Significant!
# Use Dunn post hoc to identify where significance lies
# Bonferroni's correction to account for multiple pairwise comparison 
dunnTest(dist_per_sec ~ period, data = spgr_td_2, method = "bonferroni")

# Plot distance per period - ensure it aligns with post-hoc
spgr_td_2 %>%
  filter(!is.na(period)) %>%
  ggplot() +
  geom_boxplot(aes(x = period, y = log(dist_per_sec)), fill = "indianred2", lwd = 0.75) +
  scale_x_discrete(labels = c("Dawn", "Day", "Dusk", "Night")) +
  xlab("Season") + ylab("Log distance travelled (m/s)") +
  theme_bw() +
  theme(axis.title = element_text(size = 12, face = "bold"),
        axis.text = element_text(size = 11),
        legend.title = element_text(size = 12, face = "bold"),
        legend.text = element_text(size = 11))

#----------------------
# HOME RANGE

# Can't analyse statistically - will have to do descriptive