#### All birds ####
library(ggplot2)
library(emmeans )
library(glmmTMB)
library(DHARMa)

Bird_data_reshape <- read.csv2("Data/Bird_data_reshape.csv", h= T, sep = ";")


Paludicolous_bird <- subset(Bird_data_reshape, Bird_data_reshape$Species == "Acrocephalus scirpaceus" 
                            | Bird_data_reshape$Species == "Acrocephalus arundinaceus")

Grass_bird <- subset(Bird_data_reshape, Bird_data_reshape$Species == "Alauda arvensis" | Bird_data_reshape$Species == "Galerida cristata"
                     | Bird_data_reshape$Species == "Emberiza calandra")
Edge_bird <- subset(Bird_data_reshape, Bird_data_reshape$Species == "Chloris chloris"
                    | Bird_data_reshape$Species == "Corvus corone"
                    | Bird_data_reshape$Species == "Hippolais polyglotta"
                    | Bird_data_reshape$Species == "Luscinia megarhynchos" 
                    | Bird_data_reshape$Species == "Picus viridis"
                    | Bird_data_reshape$Species == "Pica pica" 
                    | Bird_data_reshape$Species == "Sylvia atricapilla"
                    | Bird_data_reshape$Species == "Parus major"
                    | Bird_data_reshape$Species == "Upupa epops")

Paludicolous_bird$bird_type <- "Paludicolous"
Grass_bird$bird_type <- "Grass"
Edge_bird$bird_type <- "Edge"

# make one df
All_birds <- rbind(Paludicolous_bird,Grass_bird,Edge_bird)
head(All_birds)
hist(All_birds$Abundance)

# rename species (easier to plot)
All_birds$Species_short[All_birds$Species=="Acrocephalus scirpaceus"] <- "A. scirpaceus"
All_birds$Species_short[All_birds$Species=="Acrocephalus arundinaceus"] <- "A. arundinaceus"
All_birds$Species_short[All_birds$Species=="Alauda arvensis"] <- "A. arvensis"
All_birds$Species_short[All_birds$Species=="Galerida cristata"] <- "G. cristata"
All_birds$Species_short[All_birds$Species=="Emberiza calandra"] <- "E. calandra"
All_birds$Species_short[All_birds$Species=="Chloris chloris"] <- "C. chloris"
All_birds$Species_short[All_birds$Species=="Corvus corone"] <- "C. corone"
All_birds$Species_short[All_birds$Species=="Hippolais polyglotta"] <- "H. polyglotta"
All_birds$Species_short[All_birds$Species=="Luscinia megarhynchos"] <- "L. megarhynchos"
All_birds$Species_short[All_birds$Species=="Picus viridis"] <- "P. viridis"
All_birds$Species_short[All_birds$Species=="Pica pica"] <- "P. pica"
All_birds$Species_short[All_birds$Species=="Sylvia atricapilla"] <- "S. atricapilla"
All_birds$Species_short[All_birds$Species=="Parus major"] <- "P. major"
All_birds$Species_short[All_birds$Species=="Upupa epops"] <- "U. epops"

# from m2 to ha
All_birds$S_hedgerow_ha <- All_birds$S_hedgerow*0.0001
All_birds$S_grass_strip_ha <- All_birds$S_grass_strip*0.0001
All_birds$S_reeds_ha <- All_birds$S_reeds*0.0001
All_birds$S_forest_ha <- All_birds$S_forest*0.0001
All_birds$S_open_ha <- All_birds$S_open*0.0001
All_birds$S_wetland_ha <- All_birds$S_wetland*0.0001

# model
mod <- glmmTMB(Abundance ~ Species_short*scale(S_hedgerow_ha) + Species_short*scale(S_grass_strip_ha) + Species_short*scale(S_reeds_ha) 
               + Species_short*scale(S_forest_ha) + Species_short*scale(S_open_ha) + Species_short*scale(S_wetland_ha)
               + scale(Crop_MFS)*Species_short + scale(Crop_SHDI)*Species_short
               + (1|Crop)+(1|Field_number)
               + exp(pos.round + 0|group)
               ,family=nbinom2()
               ,data=All_birds)

summary(mod)
plot(allEffects(mod))
plot(simulateResiduals(mod)) 
performance::check_overdispersion(mod) 

# extract slope estimates
vec_var <- c("S_hedgerow_ha","S_grass_strip_ha","S_reeds_ha","S_forest_ha","S_open_ha","S_wetland_ha","Crop_SHDI","Crop_MFS") 
L_var <- list()
for(i in 1:length(vec_var)){
  temp <- as.data.frame(emtrends(mod, pairwise ~ Species_short, var = vec_var[i] )$emtrends)
  colnames(temp) <- c("Species_short","estimate","SE","df" , "lower" , "upper")
  temp$var <- vec_var[i]
  L_var[[i]] <- temp
}

M_var <- do.call(rbind,L_var)
head(M_var)

# add habitat preference
All_birds_dupli <- All_birds[! duplicated(All_birds$Species_short ), ]
M_var_type <- merge(M_var,All_birds_dupli[,c("Species_short","bird_type")],by="Species_short")
M_var_type$bird_type[M_var_type$bird_type=="Paludicolous"] <- "Reedbed birds"
M_var_type$bird_type[M_var_type$bird_type=="Grass"] <- "Grassland birds"
M_var_type$bird_type[M_var_type$bird_type=="Edge"] <- "Forest edge birds"

# clean habitat names
M_var_type$var[M_var_type$var=="S_hedgerow_ha"] <- "Hedgerows"
M_var_type$var[M_var_type$var=="S_grass_strip_ha"] <- "Grass strips"
M_var_type$var[M_var_type$var=="S_reeds_ha"] <- "Reed strips"
M_var_type$var[M_var_type$var=="S_forest_ha"] <- "Woodland"
M_var_type$var[M_var_type$var=="S_open_ha"] <- "Grassland"
M_var_type$var[M_var_type$var=="S_wetland_ha"] <- "Wetland"


# plot
ggplot(subset(M_var_type,M_var_type$var!="Crop_SHDI"&M_var_type$var!="Crop_MFS"))+
  geom_point(aes(x=factor(Species_short,c("C. chloris","C. corone","H. polyglotta","L. megarhynchos","P. major","P. pica","P. viridis","S. atricapilla","U. epops","A. arvensis","E. calandra","G. cristata","A. arundinaceus","A. scirpaceus")),y=estimate,color=bird_type), position=position_dodge(width=0.1), alpha = 0.7)+
  geom_errorbar(aes(x=factor(Species_short,c("C. chloris","C. corone","H. polyglotta","L. megarhynchos","P. major","P. pica","P. viridis","S. atricapilla","U. epops","A. arvensis","E. calandra","G. cristata","A. arundinaceus","A. scirpaceus")),ymin=lower,ymax=upper,col=bird_type),width=0.0001, position=position_dodge(width=0.1), alpha = 0.7)+
  scale_color_manual(values= c(pal_d3("category10")(3)[c(3,2,1)]),name="Habitat \npreference")+
  facet_wrap(~factor(var,c("Woodland","Hedgerows","Grassland","Grass strips","Wetland","Reed strips")),scales="free",nc=2)+
  ylab("Slope estimate")+
  xlab("Species")+
  geom_hline(yintercept =0)+
  theme_bw()+
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1,face="italic")) 

ggsave("plot_all_species.png", width = 8, height = 10,dpi=300) 


ggplot(subset(M_var_type,M_var_type$var=="Crop_SHDI"|M_var_type$var=="Crop_MFS"))+
  geom_point(aes(x=factor(Species_short,c("C. chloris","C. corone","H. polyglotta","L. megarhynchos","P. major","P. pica","P. viridis","S. atricapilla","U. epops","A. arvensis","E. calandra","G. cristata","A. arundinaceus","A. scirpaceus")),y=estimate,color=bird_type), position=position_dodge(width=0.1), alpha = 0.7)+
  geom_errorbar(aes(x=factor(Species_short,c("C. chloris","C. corone","H. polyglotta","L. megarhynchos","P. major","P. pica","P. viridis","S. atricapilla","U. epops","A. arvensis","E. calandra","G. cristata","A. arundinaceus","A. scirpaceus")),ymin=lower,ymax=upper,col=bird_type),width=0.0001, position=position_dodge(width=0.1), alpha = 0.7)+
  scale_color_manual(values= c(pal_d3("category10")(3)[c(3,2,1)]),name="Habitat \npreference")+
  facet_wrap(~factor(var,c("Crop_SHDI","Crop_MFS")),scales="free",nc=2)+
  ylab("Slope estimate")+
  xlab("Species")+
  geom_hline(yintercept =0)+
  theme_bw()+
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1,face="italic"))