library('agricolae')
library('tidyverse')
library('reshape2')
library('vegan')
library('ape')
library('plyr')
library('ggrepel')
library('ggthemes')
library('corrplot')
library('mgcv')
library('lmerTest')
library('MuMIn')
library('DHARMa')
library('fBasics')
library('multcomp')
library('car')
library('emmeans')
library('ggeffects')
library('glmmTMB')
library('ggplot2')
library('gridExtra')
library('Matrix')
library('writexl')
library('modEvA')
library('Rmisc')

library(reshape2)
library(dplyr)
library(ggplot2)
library(agricolae)
library(corrplot)
library(ggeffects)
library(gridExtra)
library(emmeans)
library(glmmTMB)
library(MuMIn)
library(ggthemes)


# Behavioural adaptations of the Common Buzzard Buteo buteo
# for foraging along the expressways

# Dataset loading
library(readxl)
broad <- read_xlsx('broad.xlsx') 


###############
### MODEL 1 ###
###############

# VIF checking
car::vif (glm(Habitat~TimeHigh+TimeMedium+TimeLow,
              data=broad1), family = binomial(link='logit')) # VIF < 7.2

car::vif (glm(Habitat~TimeMedium+TimeLow,
              data=broad1), family = binomial(link='logit')) # VIF < 1.1

# Model development
mod1 <- glmer(Habitat~TimeMedium+TimeLow+(1|Season),
                 data=broad1, 
                 family=binomial (link='logit'),
                 na.action=na.fail)

summary(mod1)

r.squaredGLMM (mod1)

# model prediction
library(ggeffects)
mod1_pred <- ggpredict (mod1) # model prediction

# figure 2
library(ggplot2)

Fig2 <- ggplot(data=mod1_pred[[1]], aes(x,predicted))+
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), fill="grey",alpha=0.5)+
  geom_line(linewidth=1, color="#000000")+
  geom_count(data=broad1, aes(y=Habitat, x=TimeMedium), alpha=0.5)+
  theme_classic()+
  theme (axis.title.y = element_text(size=24,vjust = +3),
         axis.title.x = element_text(size=24,vjust = -1),
         axis.text.x = element_text(size=24, color = "#000000"),
         axis.text.y = element_text(size=20, color = "#000000"),
         axis.line = element_line(linewidth=1),
         axis.ticks = element_line(linewidth = 1,color = "#000000"),
         axis.ticks.length = unit(.2, "cm"),
         legend.title = element_blank(),
         legend.text=element_text(size=20),
         legend.position = c(0.9, 0.3),
         legend.spacing.x = unit(0.3, 'cm'),
         plot.margin=unit(c(t=0.25, r=0.5, b=0.5, l=0.75), "cm"))+
  labs(y="Probability of occurrence",
       x= "Time spent on a medium-height site [s]")+
  scale_x_continuous(breaks = scales::breaks_width(100))

# figure 2 export settings
tiff(filename = "Figure2.tif", 
     unit = "cm", width=20, height=15,res = 1200, compres = "lzw")
Fig2
dev.off()

###############
### MODEL 2 ###
###############
car::vif (glm(Change~Temp+Precip+Wind,
              data=broad1), family = poisson) # VIF < 1.1

# Model development
library(lme4)

broad1$Habitat <- as.factor(broad1$Habitat)
broad1$Snow <- as.factor(broad1$Snow)
broad3$id <- as.factor(broad3$id)

library(scales)
broad1$resTemp <- rescale(broad1$Temp) 

mod2 <- glmer(Change~Habitat+Snow+Temp+Precip+Wind+(1|Season),
                  data=bdata, 
                  family=poisson,
                  na.action=na.fail) 

summary(mod2)
summary(broad1)

# Candidate models set
library(MuMIn)
mod2_dredge <- dredge(mod2) # !!! warnings

head(mod2_dredge)

# Models set within 95% confidence set limit
mod2_Avg_Mod.95 <-model.avg(mod2_dredge, cumsum(weight) <= .95)

summary(mod2_Avg_Mod.95) # <6 unique rhs>

mod2_Avg_Mod.95

# Models set that fall in 95% confidence set limit
head(mod2_Avg_Mod.95,6)


# Best model (according to AICc)
mod2_fin <- glmer(Change~Habitat+Temp+Wind+(1|Season),
              data=broad1, 
              family=poisson,
              na.action=na.fail) # !!! nie może być funkcji logit

summary(mod2_fin)

r.squaredGLMM (mod2_fin)


###############
### MODEL 3 ###
###############
car::vif (glm(Attack~Change+Temp+Precip+Wind,
              data=broad1), family = poisson) # VIF < 1.1

# Model development
mod3 <- glmer(Attack~Habitat+Snow+Temp+Precip+Wind+(1|Season),
            data=bdata1, 
            family=poisson,
            na.action=na.fail) # !!! nie może być funkcji logit

summary(mod3)

# Candidate models set
mod3_dredge <- dredge(mod3) ad(mod3_dredge)


# Models set within 95% confidence set limit
mod3_Avg_Mod.95 <-model.avg(mod3_dredge, cumsum(weight) <= .95)

summary(mod3_Avg_Mod.95) # <13 unique rhs>

mod3_Avg_Mod.95

# Models set that fall in 95% confidence set limit
head(mod3_Avg_Mod.95,13)

# Best model (according to AICc)
mod3_fin <- glmer(Attack~Habitat+Snow+Wind+(1|Season),
                  data=bdata1, 
                  family=poisson,
                  na.action=na.fail) 

summary(mod3_fin)

r.squaredGLMM (mod3_fin)

#PLOT ATTACK
library(reshape2)

summary(broad1)
colnames(broad1)

broad_att <- melt(broad1, id.vars='Habitat',  
               measure.vars=c('Attack', 'Change'))

library(Rmisc)
broad_att_Stat <- summarySE(broad_att, measurevar="value", groupvars=c("Habitat", "variable"))

library(ggplot2)
tiff(filename = "Praca3_Fig2.tif", 
     unit = "cm", width=20, height=15,res = 1200, compres = "lzw")

ggplot(broad_att_Stat, aes(x=Habitat, y=value, color=variable))+
  geom_errorbar(aes(ymin = value-ci, ymax = value+ci), width=0.2, linewidth=1, position = position_dodge(width = 0.3))+
  geom_point (aes(shape=variable), size=6, position = position_dodge(width = 0.3))+
  scale_shape_manual("", values=c(15,16), labels=c("Attack", "Change"))+
  scale_color_manual("", values=c('black','grey60'), labels=c("Attack", "Change"))+
  theme_classic()+
  theme(axis.title.x = element_text(size = 24,vjust = 0),
        axis.title.y = element_text(size = 24,vjust = +3),
        axis.text.x = element_text(size = 24, color = "#000000"),
        axis.text.y = element_text(size = 20, color = "#000000"),
        axis.ticks.y = element_line(linewidth=1),
        axis.ticks.length=unit(0.2,"cm"),
        axis.ticks.x = element_blank(),
        axis.ticks = element_line(linewidth = 0.5,color = "#000000"),
        axis.line=element_line(linewidth=1),
        legend.text=element_text(size=20),
        legend.position = c(0.2, 0.9),
        legend.spacing.x = unit(0.3, 'cm'),
        panel.border = element_blank(),
        plot.margin=unit(c(0.5,0.5,0.25,0.75), "cm"))+
  scale_y_continuous(expand = c(0, 0),limits = c(0,1.4),  breaks = seq(0,1.4, by = 0.2))+
  scale_x_discrete(labels=c("0" = "Farmland", "1" = "Expressway"))+
  labs(x = "Habitat type",
       y="Number of attacks or changes")
 
dev.off()

#PLOT SNOW
summary(broad)

broad_snow <- melt(broad1, id.vars='Snow',  
                  measure.vars=c("Attack", "Change"))

broad_snow_Stat <- summarySE(broad_snow, measurevar="value", groupvars=c("Snow","variable"))

broad_snow_Stat$snow <- as.factor(broad_snow_Stat$snow)

tiff(filename = "Praca3_Fig3.tif", 
     unit = "cm", width=20, height=15,res = 1200, compres = "lzw")

ggplot(broad_snow_Stat, aes(x=Snow, y=value, color=variable))+
  geom_errorbar(aes(ymin = value-ci, ymax = value+ci), width=0.2, linewidth=1, position = position_dodge(width = 0.3))+
  geom_point (aes(shape=variable), size=6, position = position_dodge(width = 0.3))+
  scale_shape_manual("", values=c(15,16), labels=c("Attack", "Change"))+
  scale_color_manual("", values=c('black','grey60'), labels=c("Attack", "Change"))+
  theme_classic()+
  theme(axis.title.x = element_text(size = 24,vjust = 0),
        axis.title.y = element_text(size = 24,vjust = +3),
        axis.text.x = element_text(size = 24, color = "#000000"),
        axis.text.y = element_text(size = 20, color = "#000000"),
        axis.ticks.y = element_line(linewidth=1),
        axis.ticks.length=unit(0.2,"cm"),
        axis.ticks.x = element_blank(),
        axis.ticks = element_line(linewidth = 0.5,color = "#000000"),
        axis.line=element_line(linewidth=1),
        legend.text=element_text(size=20),
        legend.position = c(0.85, 0.9),
        legend.spacing.x = unit(0.3, 'cm'),
        panel.border = element_blank(),
        plot.margin=unit(c(0.5,0.5,0.25,0.75), "cm"))+
  scale_y_continuous(expand = c(0, 0),limits = c(0,1.4),  breaks = seq(0,1.4, by = 0.2))+
  scale_x_discrete(labels=c("0" = "Absent", "1" = "Present"))+
  labs(x = "Snow cover",
       y="Number of attacks or changes")

dev.off()