rm(list = ls())
library(ggplot2)
library(lubridate)
library(scales)
library(gridExtra)
library(car)
library(MuMIn)
library(AICcmodavg)
library(lubridate)
library(ggpubr)
library(tidyverse)
library(plyr)
library(dplyr)
library(MASS)
library(emmeans)
library(kimisc)
library(texreg)
library(Hmisc)
library(magrittr)
library(ggeffects)
library(sjPlot)
library(jtools)
mpala <- read.csv("PeerJ_data.csv", header = T) 
str(mpala)
summary(mpala)
######################################################################################

#Statistical summaries
s1 <- ddply(mpala, c("Sex"), summarise,
            N = length(FGM_.ng.g._wet_weight),
            mean = mean(FGM_.ng.g._wet_weight),
            sd = sd(FGM_.ng.g._wet_weight),
            se = sd/sqrt(N))
s1

s2 <- ddply(mpala, c("Group_type"), summarise,
            N = length(FGM_.ng.g._wet_weight),
            mean = mean(FGM_.ng.g._wet_weight),
            sd = sd(FGM_.ng.g._wet_weight),
            se = sd/sqrt(N))
s2

s3 <- ddply(mpala, c("Body_condition_Score"), summarise,
            N = length(FGM_.ng.g._wet_weight),
            mean = mean(FGM_.ng.g._wet_weight),
            sd = sd(FGM_.ng.g._wet_weight),
            se = sd/sqrt(N))
s3

s4 <- ddply(mpala, c("Reaction_Index"), summarise,
            N = length(FGM_.ng.g._wet_weight),
            mean = mean(FGM_.ng.g._wet_weight),
            sd = sd(FGM_.ng.g._wet_weight),
            se = sd/sqrt(N))
s4

s5 <- ddply(mpala, c("Ranging_behavior"), summarise,
            N = length(FGM_.ng.g._wet_weight),
            mean = mean(FGM_.ng.g._wet_weight),
            sd = sd(FGM_.ng.g._wet_weight),
            se = sd/sqrt(N))
s5

s6 <- ddply(mpala, c("Age_group"), summarise,
            N = length(FGM_.ng.g._wet_weight),
            mean = mean(FGM_.ng.g._wet_weight),
            sd = sd(FGM_.ng.g._wet_weight),
            se = sd/sqrt(N))
s6

df1 <- table(mpala$Ranging_behavior, mpala$Reaction_Index)
df1
chisq.test(df1)
df2 <- cbind(df1[, "Nonresponsive within 20m"], df1[, "Retreating on vehicle approach"] + 
               df1[, "Running away"])
df2
chisq.test(df2)

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

#Fitting the basic models

m1 <- lm(log(FGM_.ng.g._wet_weight) ~ Ranging_behavior + Reaction_Index + 
           NDVI, data = mpala)
Anova(m1)
summ(m1, confint = TRUE, digits = 4)
emmeans(m1, "Reaction_Index", type = "response")
lsmeans(m1, pairwise ~ Ranging_behavior, adjust = "tukey")
lsmeans(m1, pairwise ~ Reaction_Index, adjust = "tukey")
par(mfrow = c(2, 2))
plot(m1)

#Fitting models with additional variables influencing FGM

m2 <- lm(log(FGM_.ng.g._wet_weight) ~ Ranging_behavior + Reaction_Index + 
           NDVI + Time_delay, data = mpala)
Anova(m2)
summ(m2, confint = TRUE, digits = 4)

m3 <- lm(log(FGM_.ng.g._wet_weight) ~ Ranging_behavior + Reaction_Index + 
           NDVI + Group_type, data = mpala)
Anova(m3)
summ(m3, confint = TRUE, digits = 4)
lsmeans(m3, pairwise ~ Group_type, adjust = "tukey")


m4 <- lm(log(FGM_.ng.g._wet_weight) ~ Ranging_behavior + Reaction_Index + 
           NDVI + Age_group, data = mpala)
Anova(m4)
summ(m4, confint = TRUE, digits = 4)
lsmeans(m4, pairwise ~ Age_group, adjust = "tukey")

m5 <- lm(log(FGM_.ng.g._wet_weight) ~ Ranging_behavior + Reaction_Index + 
           NDVI + Sex, data = mpala)
Anova(m5)
summ(m5, confint = TRUE, digits = 4)

m6 <- lm(log(FGM_.ng.g._wet_weight) ~ Ranging_behavior + Reaction_Index + 
           NDVI + factor(Body_condition_Score), data = mpala)
Anova(m6)
summ(m6, confint = TRUE, digits = 4)
emmeans(m6, "Body_condition_Score", type = "response")
lsmeans(m6, pairwise ~ Body_condition_Score, adjust = "tukey")


m7 <- lm(log(FGM_.ng.g._wet_weight) ~ Ranging_behavior + Reaction_Index + 
           NDVI + Total_counts, data = mpala)
Anova(m7)
summ(m7, confint = TRUE, digits = 4)

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

# Model selection based on AICc

model_list <- tibble::lst(m1, m2, m3, m4, m5, m6, m7)
aic_table <- aictab(model_list) 
aic_table


#Plotting the basic model

p1 <- ggpredict(m1, "Ranging_behavior")
p1
raw <- attr(p1, "rawdata")

gp1 <- ggplot(p1, aes(x = x, y = predicted)) + geom_point(position = position_dodge(.2)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2, position = position_dodge(.9)) + 
  ylim(3, 7) + geom_jitter(data = raw, mapping = aes(x = x, y = response, colour = "blue")) + 
  theme(legend.position = "none")
gp_1 <- gp1 + labs(x = "\nRanging behavior", y = "Predicted FGM (ng/g)", title = NULL)
gp_1 + theme_classic(base_size = 12) + theme(axis.text.x = element_text(angle = 0, size = 12, vjust = 0.4)) + 
  theme(axis.text.x = element_text(colour = "black", vjust = 0.35), axis.text.y = 
          element_text(colour = "black", vjust = 0.35)) + theme(axis.title.y = element_text(
            margin = margin(t = 0, r = 20, b = 0, l = 0)
          )) + theme(legend.position = "none") 


p2 <- ggpredict(m1, "Reaction_Index")
p2
raw2 <- attr(p2, "rawdata")
gp2 <- ggplot(p2, aes(x = x, y = predicted)) + geom_point(position = position_dodge(.2)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2, position = position_dodge(.9)) + ylim(3, 7) +
  geom_jitter(data = raw2, mapping = aes(x = x, y = response, colour = "blue")) + theme(legend.position = "none")
gp_2 <- gp2 + labs(x = "\nReaction index", y = "Predicted FGM (ng/g)", title = NULL)
gp_2 + theme_classic(base_size = 12) + theme(axis.text.x = element_text(angle = 0, size = 12, vjust = 0.4)) + 
  theme(axis.text.x = element_text(colour = "black", vjust = 0.35), axis.text.y = 
          element_text(colour = "black", vjust = 0.35)) + theme(axis.title.y = element_text(
            margin = margin(t = 0, r = 20, b = 0, l = 0)
          )) + theme(legend.position = "none")

p3 <- ggpredict(m1, "NDVI")
p3
raw3 <- attr(p3, "rawdata")
gp3 <- ggplot(raw3, aes(x = x, y = response)) + geom_point(color = "grey") + 
  geom_smooth(aes(x = x, y = response), method = "lm", formula = y ~ x + I(x^2))
gp_3 <- gp3 + labs(x = "\nNDVI", y = "\nPredicted FGM (ng/g)", title = NULL) + theme_bw() + theme_classic()
gp_3 + theme_classic(base_size = 12) + theme(axis.text.x = element_text(angle = 0, size = 12, vjust = 0.5)) + 
  theme(axis.text.x = element_text(colour = "black", vjust = 0.5), axis.text.y = 
          element_text(colour = "black", vjust = 0.5)) + theme(axis.title.y = element_text(
            margin = margin(t = 0, r = 20, b = 0, l = 0)
          ))

#################################################################################################
#Graph summary plot
Date2 <- as.Date(mpala$Date, format = "%d-%b-%y") 
class(Date2)
head(Date2) 
s6 <- ggplot(mpala, aes(x = Date2, y = NDVI)) + geom_point(color = "grey") + 
  geom_smooth(color = "black", se = FALSE, method = "gam")
s6_1 <- s6 + labs(x = "\nMonths of data collection (2019)", y = "NDVI", title = NULL) + 
  theme_bw() + theme_classic()
s6_1 + theme_classic(base_size = 12) + theme(axis.text.x = element_text(angle = 0, size = 12, vjust = 0.5)) + 
  theme(axis.text.x = element_text(colour = "black", vjust = 0.5), axis.text.y = 
          element_text(colour = "black", vjust = 0.5)) + theme(axis.title.y = element_text(
            margin = margin(t = 0, r = 20, b = 0, l = 0)
          ))

# Creating a table from the model. 

tab_model(m1, m3, m4,
          pred.labels = c("Intercept", "Ranging_behavior", "Reaction_Index", "NDVI",
                          "Group_type", "Age_group"),
          dv.labels = c("Model 1", "Model 3", "Model 4"),
          string.pred = "Coeffcient",
          string.ci = "Conf. Int (95%)",
          string.p = "P-Value", digits.p = 2)
################################################################################################