###############################################################################
#### STATISTICAL ANALYSES PART 2 - FINE-SCALE TRACKING OF GREEN-UP -- ISSA ####
###############################################################################

###-- PACKAGES --###
library(ggplot2)
library(lubridate)
library(terra)
library(amt)
library(survival)
library(raster)
library(tidyverse)
library(gsubfn)
library(doParallel)
library(sf)
library(ggcorrplot)
library(patchwork)
library(scales)
library(glmmTMB)
###--------------###

###-- DATA --###
load("E:/ibex_migration_thesis/MS3-proba_migration/MS/Ecology_letters/data_zenodo/data/data_issa.RData")
load("E:/ibex_migration_thesis/MS3-proba_migration/MS/Ecology_letters/data_zenodo/data/data_spring.RData")
###----------###

### Load data_issa.Rdata

## Each row of the dataframe corresponds to a step (two consecutive locations) from a track of an animal
## for each step, the following variables are given:
# - case: either 0 or 1, 0 if the step is a random step, 1 if the step is a real step performed by the animal
# - step_id2: step id, each real step shares a step id with its 15 associated random steps
# - id: individual id for the individual-year
# - ani_id: animal id 
# - pop: population of monitoring
# - t1_: start date of the steps
# - t2_: end date of the steps
# - dt_: timestep (hours)
# - sl_: length of the step (meters)
# - log_sl_: log(length of the step)
# - ta_: turning angle of the step
# - cos_ta_: cosinus(turning angle)
# - DFP: Days-from-peak value at the end of the step
# - IRG: IRG value at the end of the step
# - north: Northness value (between -1 and 1, 1 for north) at the end of the step
# - slope: Slope value at the end of the step 

##---------------------------------------------------------------------------------


### Load data_spring.Rdata

## Each row of the dataframe corresponds to an individual-year 
## for each individual-year the following variables are given:
# - CIRGs: average IRG at locations of ibex on their summer range (cf Methods, section 5: Exposure to green-up)
# - DFP: average Days-from-peak at locations of ibex on their summer range (cf Methods, section 5: Exposure to green-up)
# - TIRG: Time exposed to IRG values >0.4 (number of days) on the summer range of ibex (cf Methods, section 5: Exposure to green-up)
# - SpeedGU: Values of speed of green-up (depends on the year and popualtion of monitoring) (cf Methods, section 4: Spring vegetation dynamics and springness metrics)
# - SHGU: Values of spatial heterogeneity of green-up (depends on the year and popualtion of monitoring) (cf Methods, section 4: Spring vegetation dynamics and springness metrics)
# - SpeedGU_WP: Temporal component of the speed of green-up variable (SpeedGU - population mean) (cf Methods, section 4: Spring vegetation dynamics and springness metrics)
# - SpeedGU_BP: Spatial component of the speed of green-up variable (population mean) (cf Methods, section 4: Spring vegetation dynamics and springness metrics)
# - SHGU_WP: Temporal component of the spatial heterogeneity of green-up variable (SHGU - population mean) (cf Methods, section 4: Spring vegetation dynamics and springness metrics)
# - SHGU_BP: Spatial component of the spatial heterogeneity of green-up variable (population mean) (cf Methods, section 4: Spring vegetation dynamics and springness metrics)
# - cat1: Category of speed of green-up based on the quantiles 0.33 and 0.66 of SpeedGU_WP, either low, average or high
# - cat2: Category of spatial heterogeneity of green-up based on the quantiles 0.33 and 0.66 of SHGU_WP, either low, average or high
# - year: year of monitoring
# - pop: population of monitoring
# - id: individual id for the individual-year
# - new_id: animal id 
# - bhvr1: migratory behaviour - either migrant or resident
# - AgeCl: Age class either 1, 2 or 3 (cf Methods, section 2: Individual characteristics)
# - sex: sex of the animal 



## The following analyses consist of integrated Step Selection Analyses (conditional logistic regression) performed 
## to assess the fine scale selection for spring green-up (IRG or Days-from-peak)

## ONE MODEL WILL BE FITTED FOR EACH COMBINATION OF SPRING TYPE (cat1 and  cat2)/SEX/MIGRATORY BEHAVIOUR
table(data_spring$cat1, data_spring$sex, data_spring$bhvr1)
table(data_spring$cat2, data_spring$sex, data_spring$bhvr1)


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

### 1) FIT MODELS FOR THE DIFFERENT CATEGORIES OF SPEED OF GREEN UP (cat1)

OUT2_issf_cat1 <- NULL
for(c in unique(data_spring$cat1)){
  print(c)
  for(s in c("M", "F")){
    print(s)
    for(b in c("migrant", "resident")){
      print(b)
      id_cat <- data_spring$id[data_spring$cat1 %in% c & data_spring$sex %in% s & data_spring$bhvr1 %in% b]
      nb <- length(id_cat)
      
      dfk_c <- data_issa[data_issa$id %in% id_cat,]
      
      #removing steps in forested areas
      sid <- unique(dfk_c$step_id2[dfk_c$foret %in% 1 & dfk_c$case_ %in% "TRUE"])
      dfk_c <- dfk_c[!dfk_c$step_id2 %in% sid,]
      
      moddata <- dfk_c %>% dplyr::select(c(case, DFP, IRG, north, slope, sl_, log_sl_, cos_ta_, step_id2, ani_id, pop))
      moddata <- na.omit(moddata)
      
      ## Fit ISSA with IRG as the environmental covariate
      TMB_IRG <- glmmTMB(case ~ -1 + log_sl_ + sl_ + cos_ta_ + north + slope +  IRG +(1|step_id2) + 
                           (1|ani_id) + (1|pop),
                         family=poisson, data = moddata,
                         start = list(theta=c(log(1e3), rep(0,2))),  
                         map = list(theta=factor(c(NA, 1:2))) )
      
      ## Fit ISSA with DFP as the environmental covariate
      TMB_DFP <- glmmTMB(case ~ -1 + log_sl_ + sl_ + cos_ta_ + north + slope +  DFP +(1|step_id2) + 
                           (1|ani_id) + (1|pop),
                         family=poisson, data = moddata,
                         start = list(theta=c(log(1e3), rep(0,2))),
                         map = list(theta=factor(c(NA, 1:2))) )
      
      res1 <- data.frame(summary(TMB_IRG)$coefficients$cond)
      res2 <- data.frame(summary(TMB_DFP)$coefficients$cond)
      colnames(res1)<-c("Coef" , "SE", "z value"  ,  "P_val" )
      colnames(res2)<-c("Coef" , "SE", "z value"  ,  "P_val" )
      res1$var<-row.names(res1)
      res1$mod <- "irg"
      res1$sex <- s
      res1$cat <- c
      res1$bhvr <- b
      res1$nb <- nb
      
      res2$var<-row.names(res2)
      res2$mod <- "dfp"
      res2$sex <- s
      res2$cat <- c
      res2$bhvr <- b
      res2$nb <- nb
      
      OUT2_issf_cat1 <- rbind(OUT2_issf_cat1, res1, res2)
      
    }
  }
}

OUT2_issf_cat1$SS[OUT2_issf_cat1$P_val < 0.05 & OUT2_issf_cat1$Coef < 0] <- "Neg"
OUT2_issf_cat1$SS[OUT2_issf_cat1$P_val < 0.05 & OUT2_issf_cat1$Coef > 0] <- "Pos"
OUT2_issf_cat1$SS[OUT2_issf_cat1$P_val >= 0.05] <- "NS"

OUT2_issf_cat1 <- OUT2_issf_cat1[OUT2_issf_cat1$var %in% c("IRG", "DFP"),]

## Compute Confidence Intervals
OUT2_issf_cat1$ICup <- OUT2_issf_cat1$Coef + 1.96*OUT2_issf_cat1$SE
OUT2_issf_cat1$ICinf <- OUT2_issf_cat1$Coef - 1.96*OUT2_issf_cat1$SE

## Order factor levels for plotting results
OUT2_issf_cat1$cat <- ordered(OUT2_issf_cat1$cat, levels = c("Low", "Average", "High"))
OUT2_issf_cat1$sex <- as.factor(OUT2_issf_cat1$sex)
OUT2_issf_cat1$sex <- relevel(OUT2_issf_cat1$sex, ref = c("M"))

##------ FIGURES 4A and 4C IN MANUSCRIPT ------##
## IRG - FIG 4A
gc1 <-ggplot(OUT2_issf_cat1[OUT2_issf_cat1$var %in% "IRG",], aes(x=cat, y=Coef, col=bhvr)) + 
  geom_hline(aes(yintercept=0), col="black", linetype=2)+
  geom_point(position=position_dodge(width=0.5), size=5) + facet_grid(. ~ sex) +
  geom_linerange(aes(ymin=ICinf, ymax=ICup, col=bhvr), linewidth=1, position=position_dodge(width=0.5)) + 
  theme_classic() + ylim(-0.2,0.8) +
  theme(axis.text=element_text(size=26),
        strip.text= element_text(size = 20, face="bold"),
        axis.text.x = element_text(angle=0, vjust=0.8, hjust=0.8, size=20),
        axis.title=element_text(size=26),
        legend.title = element_text(size = 24, face='bold'),
        legend.text = element_text(size=24),
        #legend.position = c(2,2),
        legend.position = "none",
        plot.margin = margin(10, 30, 10, 10),
        text = element_text(family = "serif"),
        axis.ticks.length=unit(.3, "cm"),
        strip.text.x = element_blank(),
        panel.spacing = unit(2, "cm", data = NULL)) +
  xlab("") +
  ylab("Selection coefficient IRG") +
  scale_color_manual(values=c("firebrick","dodgerblue4"))+
  guides(fill=guide_legend(""), color =guide_legend("") , shape= guide_legend(""), linetype=guide_legend(""))

## DFP - FIG 4C
gc2 <-ggplot(OUT2_issf_cat1[OUT2_issf_cat1$var %in% "DFP",], aes(x=cat, y=Coef, col=bhvr)) + 
  geom_hline(aes(yintercept=0), col="black", linetype=2)+
  geom_point(position=position_dodge(width=0.5), size=5) + facet_grid(. ~ sex) +
  geom_linerange(aes(ymin=ICinf, ymax=ICup, col=bhvr), linewidth=1, position=position_dodge(width=0.5)) + 
  theme_classic() + ylim(-0.03,0) +
  theme(axis.text=element_text(size=26),
        strip.text= element_text(size = 20, face="bold"),
        axis.text.x = element_text(angle=0, vjust=0.8, hjust=0.8, size=20),
        axis.title=element_text(size=26),
        legend.title = element_text(size = 24, face='bold'),
        legend.text = element_text(size=24),
        #legend.position = c(2,2),
        legend.position = "none",
        plot.margin = margin(10, 30, 10, 10),
        text = element_text(family = "serif"),
        axis.ticks.length=unit(.3, "cm"),
        strip.text.x = element_blank(),
        panel.spacing = unit(2, "cm", data = NULL)) +
  xlab("Speed of green-up") +
  ylab("Selection coefficient DFP") +
  scale_color_manual(values=c("firebrick","dodgerblue4"))+
  guides(fill=guide_legend(""), color =guide_legend("") , shape= guide_legend(""), linetype=guide_legend(""))


gc1/gc2


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

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

### 2) FIT MODELS FOR THE DIFFERENT CATEGORIES OF SPATIAL HETEROGENEITY OF GREEN-UP (cat2)


OUT2_issf_cat2 <- NULL
for(c in unique(data_spring$cat2)){
  print(c)
  for(s in c("M", "F")){
    print(s)
    for(b in c("migrant", "resident")){
      print(b)
      id_cat <- data_spring$id[data_spring$cat2 %in% c & data_spring$sex %in% s & data_spring$bhvr1 %in% b]
      nb <- length(id_cat)
      
      dfk_c <- data_issa[data_issa$id %in% id_cat,]
      
      #removing steps in forested areas
      sid <- unique(dfk_c$step_id2[dfk_c$foret %in% 1 & dfk_c$case_ %in% "TRUE"])
      dfk_c <- dfk_c[!dfk_c$step_id2 %in% sid,]
      
      moddata <- dfk_c %>% dplyr::select(c(case, DFP, IRG, north, slope, sl_, log_sl_, cos_ta_, step_id2, ani_id, pop))
      moddata <- na.omit(moddata)
      
      moddata <- moddata %>% group_by(step_id2) %>% mutate(verif = length(unique(moddata$case[moddata$step_id2 %in% step_id2])))
      
      ## Fit ISSA with IRG as the environmental covariate
      TMB_IRG <- glmmTMB(case ~ -1 + log_sl_ + sl_ + cos_ta_ + north + slope + IRG +(1|step_id2) +
                           (1|ani_id) + (1|pop),
                         family=poisson, data = moddata,
                         start = list(theta=c(log(1e3), rep(0,2))),  
                         map = list(theta=factor(c(NA, 1:2))) )
      
      
      ## Fit ISSA with DFP as the environmental covariate
      TMB_DFP <- glmmTMB(case ~ -1 + log_sl_ + sl_ + cos_ta_ + north + slope + DFP +(1|step_id2) +
                           (1|ani_id) + (1|pop),
                         family=poisson, data = moddata,
                         start = list(theta=c(log(1e3), rep(0,2))),   
                         map = list(theta=factor(c(NA, 1:2))) )
      
      res1 <- data.frame(summary(TMB_IRG)$coefficients$cond)
      res2 <- data.frame(summary(TMB_DFP)$coefficients$cond)
      colnames(res1)<-c("Coef" , "SE", "z value"  ,  "P_val" )
      colnames(res2)<-c("Coef" , "SE", "z value"  ,  "P_val" )
      res1$var<-row.names(res1)
      res1$mod <- "irg"
      res1$sex <- s
      res1$cat <- c
      res1$bhvr <- b
      res1$nb <- nb
      
      res2$var<-row.names(res2)
      res2$mod <- "dfp"
      res2$sex <- s
      res2$cat <- c
      res2$bhvr <- b
      res2$nb <- nb
      
      OUT2_issf_cat2 <- rbind(OUT2_issf_cat2, res1, res2)
      
    }
  }
}


OUT2_issf_cat2$SS[OUT2_issf_cat2$P_val < 0.05 & OUT2_issf_cat2$Coef < 0] <- "Neg"
OUT2_issf_cat2$SS[OUT2_issf_cat2$P_val < 0.05 & OUT2_issf_cat2$Coef > 0] <- "Pos"
OUT2_issf_cat2$SS[OUT2_issf_cat2$P_val >= 0.05] <- "NS"

OUT2_issf_cat2 <- OUT2_issf_cat2[OUT2_issf_cat2$var %in% c("IRG", "DFP"),]

## Compute Confidence Intervals
OUT2_issf_cat2$ICup <- OUT2_issf_cat2$Coef + 1.96*OUT2_issf_cat2$SE
OUT2_issf_cat2$ICinf <- OUT2_issf_cat2$Coef - 1.96*OUT2_issf_cat2$SE

## Order factor levels for plotting results
OUT2_issf_cat2$cat <- ordered(OUT2_issf_cat2$cat, levels = c("Low", "Average", "High"))
OUT2_issf_cat2$sex <- as.factor(OUT2_issf_cat2$sex)
OUT2_issf_cat2$sex <- relevel(OUT2_issf_cat2$sex, ref = c("M"))


##------ FIGURES 4B and 4D IN MANUSCRIPT ------##
## IRG - FIG 4B
gc3 <-ggplot(OUT2_issf_cat2[OUT2_issf_cat2$var %in% "IRG",], aes(x=cat, y=Coef, col=bhvr)) + 
  geom_hline(aes(yintercept=0), col="black", linetype=2)+
  geom_point(position=position_dodge(width=0.5), size=5) + facet_grid(. ~ sex) +
  geom_linerange(aes(ymin=ICinf, ymax=ICup, col=bhvr), linewidth=1, position=position_dodge(width=0.5)) + 
  theme_classic() + ylim(-0.2,0.8) +
  theme(axis.text=element_text(size=26),
        strip.text= element_text(size = 20, face="bold"),
        axis.text.x = element_text(angle=0, vjust=0.8, hjust=0.8, size=20),
        axis.title=element_text(size=26),
        legend.title = element_text(size = 24, face='bold'),
        legend.text = element_text(size=24),
        #legend.position = c(2,2),
        legend.position = "none",
        plot.margin = margin(10, 30, 10, 10),
        text = element_text(family = "serif"),
        axis.ticks.length=unit(.3, "cm"),
        strip.text.x = element_blank(),
        panel.spacing = unit(2, "cm", data = NULL)) +
  xlab("") +
  ylab("Selection coefficient IRG") +
  scale_color_manual(values=c("firebrick","dodgerblue4"))+
  guides(fill=guide_legend(""), color =guide_legend("") , shape= guide_legend(""), linetype=guide_legend(""))


## DFP - FIG 4D
gc4 <-ggplot(OUT2_issf_cat2[OUT2_issf_cat2$var %in% "DFP",], aes(x=cat, y=Coef, col=bhvr)) + 
  geom_hline(aes(yintercept=0), col="black", linetype=2)+
  geom_point(position=position_dodge(width=0.5), size=5) + facet_grid(. ~ sex) +
  geom_linerange(aes(ymin=ICinf, ymax=ICup, col=bhvr), linewidth=1, position=position_dodge(width=0.5)) + 
  theme_classic() + ylim(-0.03,0) +
  theme(axis.text=element_text(size=26),
        strip.text= element_text(size = 20, face="bold"),
        axis.text.x = element_text(angle=0, vjust=0.8, hjust=0.8, size=20),
        axis.title=element_text(size=26),
        legend.title = element_text(size = 24, face='bold'),
        legend.text = element_text(size=24),
        legend.position = c(2,2),
        #legend.position = "none",
        plot.margin = margin(10, 30, 10, 10),
        text = element_text(family = "serif"),
        axis.ticks.length=unit(.3, "cm"),
        strip.text.x = element_blank(),
        panel.spacing = unit(2, "cm", data = NULL)) +
  xlab("Spatial-heterogeneity of green-up") +
  ylab("Selection coefficient DFP") +
  scale_color_manual(values=c("firebrick","dodgerblue4"))+
  guides(fill=guide_legend(""), color =guide_legend("") , shape= guide_legend(""), linetype=guide_legend(""))


gc3/gc4

## FIGURE 4 
(gc1/gc2)|(gc3/gc4)

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