---
title: "sliding windows and linear models"
author: "Carole Niffenegger"
date: "`r Sys.Date()`"
output: html_document
---



##### Required packages

```{r}

library(lubridate)
library(birdring)
library(ggplot2)
library(ppcor)
library(gridExtra)
library(lemon)
library(brms)
library(dplyr)

```

##### SLIDING WINDOWS

#### Correlation coefficient of start of breeding

```{r}

## CALCULATE THE  CORRELATIONS COEFFICIENTS FOR START OF BREEDING SEASON

# here, I use the normal pearson correlation coefficient rather than the partial correlation coefficient because
# I can account for the effect of the other variables in the linear model afterwards and avoid unrealistic
# combinations of variables

df.combinations.start <- expand.grid(seq(1,173, by = 2), seq(7,180, by =2))
colnames(df.combinations.start) <- c("start", "end")

# remove cases where start is after end or when time window is shorter than one week

df.combinations.start <- df.combinations.start[-which(df.combinations.start$end - df.combinations.start$start < 7),]

# create columns to fill with partial correlation coefficients

df.combinations.start$start.pcor.prec <- NA
df.combinations.start$start.pcor.prec.SE <- NA
df.combinations.start$start.pcor.temp <- NA
df.combinations.start$start.pcor.temp.SE <- NA


# load precipitation and temperature data

df.prec.year <-  read.csv("prec_year_Breg.csv", sep = ",", check.names = FALSE)
df.temp.year <- read.csv("temp_year_Breg.csv", check.names = FALSE)


df.SF <- read.csv("df.SF.csv")


# only use years between 2006 and 2021 in all datasets

df.SF <- df.SF[df.SF$year >= 2006 & df.SF$year <= 2021,]

ind.cor.start <- 1

for (ind.cor.start in 1:nrow(df.combinations.start)) {
  
  # calculate precipitation sum for one time window
  for.window <- df.combinations.start[ind.cor.start,1:2]
  # extract the doy of years of the window
  for.doy <- paste("doy", (for.window) , sep = "" )
  # extract the position of the day of year in the df.SF data file
  
  for.start.prec <- which(colnames(df.prec.year) == for.doy[1])
  for.end.prec <- which(colnames(df.prec.year) == for.doy[2])
  
  for.start.temp <- which(colnames(df.temp.year) == for.doy[1])
  for.end.temp <- which(colnames(df.temp.year) == for.doy[2])
  
  # calculate the sum of the precipitation in a given window
  
  
  prec.doy <- rowSums(df.prec.year[for.start.prec:for.end.prec]>0, na.rm = T) # number of days with precipitation
  mean.temp <- rowMeans(df.temp.year[for.start.temp:for.end.temp])

  
  for.df.mod <- data.frame(b.start = df.SF$start,
                       prec.doy = prec.doy,
                       temp.mean = mean.temp,
                       snow.timing = df.SF$snow50,
                       melt.duration = df.SF$melt.duration,
                       year = df.SF$year,
                       breg = df.SF$bioreg)
  

  ## calculate partial correlation for the length
  
  cor.prec <- cor(for.df.mod$b.start, for.df.mod$prec.doy)
  cor.temp <- cor(for.df.mod$b.start, for.df.mod$temp.mean)
  
  for.partial.SE.1 <- sqrt(1-((cor.prec)^2)) /  sqrt(nrow(for.df.mod)-2)
  for.partial.SE.2 <- sqrt(1-((cor.temp)^2)) /  sqrt(nrow(for.df.mod)-2)
  

  df.combinations.start$start.pcor.prec[ind.cor.start] <-  cor.prec
  df.combinations.start$start.pcor.temp[ind.cor.start] <-  cor.temp
  
  df.combinations.start$start.pcor.prec.SE[ind.cor.start] <- for.partial.SE.1
  df.combinations.start$start.pcor.temp.SE[ind.cor.start] <- for.partial.SE.2
  
}

df.combinations.start$start.pcor.prec.CI.low <- df.combinations.start$start.pcor.prec - (1.95 * df.combinations.start$start.pcor.prec.SE)
df.combinations.start$start.pcor.prec.CI.upr <- df.combinations.start$start.pcor.prec + (1.95 * df.combinations.start$start.pcor.prec.SE)


#write.csv(df.combinations.start, paste0("pcor_sliding_windows_start_", Sys.Date() ,".csv")) # save the output of sliding windows (including date in the file name)

df.combinations.start <- df.combinations.start[df.combinations.start$start>=91,] # reduce data to windows that start on the 1 of April or later
df.combinations.start <- df.combinations.start[df.combinations.start$end>=98,] # reduce data to windows that start on the 1 of january or later

df.combinations.start$start.date <- as.Date(origin = "2006-01-01", df.combinations.start$start)
df.combinations.start$end.date <- as.Date(origin = "2006-01-01", df.combinations.start$end)

# plot results of sliding window - precipitation

df.mask <- df.combinations.start[df.combinations.start$start.pcor.prec <= 0.25 & df.combinations.start$start.pcor.prec >= -0.25,]
min.point <- df.combinations.start[df.combinations.start$start.pcor.prec == (max(df.combinations.start$start.pcor.prec)), ]

ggplot(data = df.combinations.start, aes(x = start.date, y = end.date, col = start.pcor.prec))+
  geom_point(shape = 15, size = 3)+  
  scale_x_date(date_breaks = "2 weeks", date_labels = "%d. %b", minor_breaks = "week",  expand = c(0, 0))+
  scale_y_date(date_breaks = "2 weeks", date_labels = "%d. %b", expand = c(0, 0))+
  labs(title="Start - precipitation days")+
  scale_color_gradientn(colours = c("red", "orange", "white","lightblue", "royalblue4"), values = scales::rescale(c(-0.5,-0.25, 0, 0.25, 0.5)), limits = c(-0.7,0.7), breaks = c(-0.5,-0.25, 0, 0.25,0.5))+
  ylab("end of window")+ xlab("start of window")+
  labs(colour = "correlation coefficient", fill = "default", x = "", y = "")+
  geom_point(data = df.mask, aes(x = start.date, y = end.date), col = "grey", alpha = 0.5, shape = 15, size = 3)+
  geom_hline(yintercept = as.numeric(as.Date("2017-05-23")), linetype = "dashed")+ geom_vline(xintercept = as.numeric(as.Date("2017-05-23")), linetype = "dashed")+
  geom_point(data=min.point, aes(x = start.date, y = end.date), shape = 0, col = "black", size = 3)+
  ylab("end of window")+ xlab("start of window")+
  theme_bw()

# plot results of sliding window - temperature

df.mask <- df.combinations.start[df.combinations.start$start.pcor.temp <= 0.25 & df.combinations.start$start.pcor.temp >= -0.25,]
min.point <- df.combinations.start[df.combinations.start$start.pcor.temp == (min(df.combinations.start$start.pcor.temp)), ]

ggplot(data = df.combinations.start, aes(x = start.date, y = end.date, col = start.pcor.temp))+
  geom_point(shape = 15, size = 3)+  
  scale_x_date(date_breaks = "2 weeks", date_labels = "%d. %b", minor_breaks = "week",  expand = c(0, 0))+
  scale_y_date(date_breaks = "2 weeks", date_labels = "%d. %b", expand = c(0, 0))+
  labs(title="Start - mean temperature")+
  scale_color_gradientn(colours = c("red", "orange", "white","lightblue", "royalblue4"), values = scales::rescale(c(-0.5,-0.25, 0, 0.25, 0.5)), limits = c(-0.7,0.7), breaks = c(-0.5,-0.25, 0, 0.25,0.5))+
  ylab("end of window")+ xlab("start of window")+
  labs(colour = "correlation coefficient", fill = "default", x = "", y = "")+
  geom_point(data = df.mask, aes(x = start.date, y = end.date), col = "grey", alpha = 0.5, shape = 15, size = 3)+
  geom_hline(yintercept = as.numeric(as.Date("2017-05-23")), linetype = "dashed")+ geom_vline(xintercept = as.numeric(as.Date("2017-05-23")), linetype = "dashed")+
  geom_point(data=min.point, aes(x = start.date, y = end.date), shape = 0, col = "black", size = 3)+
  ylab("end of window")+ xlab("start of window")+
  theme(legend.position="bottom")+
  theme_bw()


```


#### Correlation coefficient of end of breeding

```{r}

## CALCULATE THE CORRELATIONS COEFFICIENT FOR END OF BREEDING SEASON

# here, I use the normal pearson correlation coefficient rather than the partial correlation coefficient because
# I can account for the effect of the other variables in the linear model afterwards and avoid unrealistic
# combinations of variables throu

df.combinations.end <- expand.grid(seq(1,220, by = 2), seq(7,227, by =2))
colnames(df.combinations.end) <- c("start", "end")

# remove cases where start is after end or when time window is shorter than one week

df.combinations.end <- df.combinations.end[-which(df.combinations.end$end - df.combinations.end$start < 7),]

# create columns to fill with partial correlation coefficients

df.combinations.end$start.pcor.prec <- NA
df.combinations.end$start.pcor.prec.SE <- NA
df.combinations.end$start.pcor.temp <- NA
df.combinations.end$start.pcor.temp.SE <- NA


# load precipitation and temperature data

df.prec.year <-  read.csv("prec_year_Breg.csv", check.names = FALSE)
df.temp.year <- read.csv("temp_year_Breg.csv", check.names = FALSE)


# only use years between 2006 and 2021 in all datasets

df.prec.year <- df.prec.year[df.prec.year$year >= 2006 & df.prec.year$year <= 2021,]
df.temp.year <- df.temp.year[df.temp.year$year >= 2006 & df.temp.year$year <= 2021,]
df.SF <- df.SF[df.SF$year >= 2006 & df.SF$year <= 2021,]



ind.cor.start <- 1 # create dummy index to test if for loop is working properly


for (ind.cor.start in 1:nrow(df.combinations.end)) {
  
  # calculate precipitation sum for one time window
  for.window <- df.combinations.end[ind.cor.start,1:2]
  # extract the doy of years of the window
  for.doy <- paste("doy", (for.window) , sep = "" )
  # extract the position of the day of year in the df.SF data file
  
  for.start.prec <- which(colnames(df.prec.year) == for.doy[1])
  for.end.prec <- which(colnames(df.prec.year) == for.doy[2])
  
  for.start.temp <- which(colnames(df.temp.year) == for.doy[1])
  for.end.temp <- which(colnames(df.temp.year) == for.doy[2])
  
  # calculate the sum of the precipitation in a given window
  
  
  prec.doy <- rowSums(df.prec.year[for.start.prec:for.end.prec]>0, na.rm = T) # number of days with precipitation
  mean.temp <- rowMeans(df.temp.year[for.start.temp:for.end.temp])

  
  for.df.mod <- data.frame(b.end = df.SF$end,
                       prec.doy = prec.doy,
                       temp.mean = mean.temp,
                       snow.timing = df.SF$snow50,
                       melt.duration = df.SF$melt.duration,
                       year = df.SF$year,
                       breg = df.SF$bioreg)
  

  ## calculate partial correlation for the length
  
  cor.prec <- cor(for.df.mod$b.end, for.df.mod$prec.doy)
  cor.temp <- cor(for.df.mod$b.end, for.df.mod$temp.mean)
  
  for.partial.SE.1 <- sqrt(1-((cor.prec)^2)) /  sqrt(nrow(for.df.mod)-2)
  for.partial.SE.2 <- sqrt(1-((cor.temp)^2)) /  sqrt(nrow(for.df.mod)-2)
  

  df.combinations.end$start.pcor.prec[ind.cor.start] <-  cor.prec
  df.combinations.end$start.pcor.temp[ind.cor.start] <-  cor.temp
  
  df.combinations.end$start.pcor.prec.SE[ind.cor.start] <- for.partial.SE.1
  df.combinations.end$start.pcor.temp.SE[ind.cor.start] <- for.partial.SE.2
  
}

df.combinations.end$start.pcor.prec.CI.low <- df.combinations.end$start.pcor.prec - (1.95 * df.combinations.end$start.pcor.prec.SE)
df.combinations.end$start.pcor.prec.CI.upr <- df.combinations.end$start.pcor.prec + (1.95 * df.combinations.end$start.pcor.prec.SE)


#write.csv(df.combinations.start, paste0("pcor_sliding_windows_start_", Sys.Date() ,".csv")) # save the output of sliding windows (including date in the file name)

df.combinations.end <- df.combinations.end[df.combinations.end$start>=121,] # reduce data to windows that start on the 15 of April or later
df.combinations.end <- df.combinations.end[df.combinations.end$end>=128,] # reduce data to windows that start on the 1 of january or later

df.combinations.end$start.date <- as.Date(origin = "2006-01-01", df.combinations.end$start)
df.combinations.end$end.date <- as.Date(origin = "2006-01-01", df.combinations.end$end)

# plot results of sliding window - precipitation

df.mask <- df.combinations.end[df.combinations.end$start.pcor.prec <= 0.25 & df.combinations.end$start.pcor.prec >= -0.25,]
min.point <- df.combinations.end[df.combinations.end$start.pcor.prec == (min(df.combinations.end$start.pcor.prec)), ]

ggplot(data = df.combinations.end, aes(x = start.date, y = end.date, col = start.pcor.prec))+
  geom_point(shape = 15, size = 2)+  
  scale_x_date(date_breaks = "2 weeks", date_labels = "%d. %b", minor_breaks = "week",  expand = c(0, 0))+
  scale_y_date(date_breaks = "2 weeks", date_labels = "%d. %b", expand = c(0, 0))+
  labs(title="End - precipitation days")+
  scale_color_gradientn(colours = c("red", "orange", "white","lightblue", "royalblue4"), values = scales::rescale(c(-0.5,-0.25, 0, 0.25, 0.5)), limits = c(-0.7,0.7), breaks = c(-0.5,-0.25, 0, 0.25,0.5))+
  ylab("end of window")+ xlab("start of window")+
  labs(colour = "correlation coefficient", fill = "default", x = "", y = "")+
  geom_point(data = df.mask, aes(x = start.date, y = end.date), col = "grey", alpha = 0.5, shape = 15, size = 2)+
  geom_hline(yintercept = as.numeric(as.Date("2017-05-23")), linetype = "dashed")+ geom_vline(xintercept = as.numeric(as.Date("2017-05-23")), linetype = "dashed")+
  ylab("end of window")+ xlab("start of window")+
  theme_bw()

# plot results of sliding window - temperature

df.mask <- df.combinations.end[df.combinations.end$start.pcor.temp <= 0.25 & df.combinations.end$start.pcor.temp >= -0.25,]
min.point <- df.combinations.end[df.combinations.end$start.pcor.temp == (min(df.combinations.end$start.pcor.temp)), ]

ggplot(data = df.combinations.end, aes(x = start.date, y = end.date, col = start.pcor.temp))+
  geom_point(shape = 15, size = 2)+  
  scale_x_date(date_breaks = "2 weeks", date_labels = "%d. %b", minor_breaks = "week",  expand = c(0, 0))+
  scale_y_date(date_breaks = "2 weeks", date_labels = "%d. %b", expand = c(0, 0))+
  labs(title="End - mean temperature")+
  scale_color_gradientn(colours = c("red", "orange", "white","lightblue", "royalblue4"), values = scales::rescale(c(-0.5,-0.25, 0, 0.25, 0.5)), limits = c(-0.7,0.7), breaks = c(-0.5,-0.25, 0, 0.25,0.5))+
  ylab("end of window")+ xlab("start of window")+
  labs(colour = "correlation coefficient", fill = "default", x = "", y = "")+
  geom_point(data = df.mask, aes(x = start.date, y = end.date), col = "grey", alpha = 0.5, shape = 15, size = 2)+
  geom_hline(yintercept = as.numeric(as.Date("2017-05-23")), linetype = "dashed")+ geom_vline(xintercept = as.numeric(as.Date("2017-05-23")), linetype = "dashed")+
  geom_point(data=min.point, aes(x = start.date, y = end.date), shape = 0, col = "black", size = 2)+
  ylab("end of window")+ xlab("start of window")+
  theme(legend.position="bottom")+
  theme_bw()


```

##### LINEAR REGRESSION (phenology and environement)


### Linear regression based on sliding windows

```{r}


### START OF THE BREEDING SEASON
### create a dataframe with evironemental variables for the time window with the max positive correlation coefficient

# time window for temperature

start.temp <- df.combinations.start$start[df.combinations.start$start.pcor.temp == min(df.combinations.start$start.pcor.temp)]
end.temp <- df.combinations.start$end[df.combinations.start$start.pcor.temp == min(df.combinations.start$start.pcor.temp)]


vec.doy.temp <- paste("doy", c(start.temp, end.temp) , sep = "" )


ind.start.temp <- which(colnames(df.temp.year) == vec.doy.temp[1])
ind.end.temp <- which(colnames(df.temp.year) == vec.doy.temp[2])

mean.temp <- rowMeans(df.temp.year[ind.start.temp:ind.end.temp])

# time window for precipitation

start.prec <- df.combinations.start$start[df.combinations.start$start.pcor.prec == max(df.combinations.start$start.pcor.prec)]
end.prec <- df.combinations.start$end[df.combinations.start$start.pcor.prec == max(df.combinations.start$start.pcor.prec)]


vec.doy.prec <- paste("doy", c(start.prec, end.prec) , sep = "" )


ind.start.prec<- which(colnames(df.prec.year) == vec.doy.prec[1])
ind.end.prec <- which(colnames(df.prec.year) == vec.doy.prec[2])

prec.doy <- rowSums(df.prec.year[ind.start.prec:ind.end.prec]>0, na.rm = T)

  
df.mod.all.windows.start <- data.frame(b.length = df.SF$length.mean,
                     b.start = df.SF$start,
                     b.end = df.SF$end,
                       b.intensity = df.SF$intensity,
                       prec.doy = prec.doy,
                       temp.mean = mean.temp,
                     snow.timing = df.SF$snow50,
                       melt.duration = df.SF$melt.duration,
                       year = df.SF$year,
                       breg = df.SF$bioreg)


mod.start <- brm(b.start ~ temp.mean * prec.doy + snow.timing +  (1|year) + (1|breg), 
                 data = df.mod.all.windows.start)
summary(mod.start)
conditional_effects(mod.start)


#saveRDS(mod.start, "mod.start.rds")
mod.start<- readRDS("mod.start.rds")


# create a plot of the marginal effect of interaction between temperature and precipitation

newdat <- expand.grid(melt.duration = mean(df.mod.all.windows.start$melt.duration),
                     snow.timing = mean(df.mod.all.windows.start$snow.timing),
                     temp.mean = seq(min(df.mod.all.windows.start$temp.mean),
                                     to = max(df.mod.all.windows.start$temp.mean),length.out =100),
                     prec.doy = c(4,11))
# draw samples from the posterior distribution based on the newly created data
pred <- posterior_epred(mod.start, newdata = newdat, re_formula = NA)


# calculate the median and uncertainty intervall from the posterior draws
newdat$pred.median <- apply(pred, 2, mean)
newdat$pred.lwr <- apply(pred, 2, quantile, 0.025)
newdat$pred.upr <- apply(pred, 2, quantile, 0.975)


newdat$pred.median.date <- as.Date(newdat$pred.median, origin = "2010-01-01")
newdat$pred.lwr.date <- as.Date(newdat$pred.lwr, origin = "2010-01-01")
newdat$pred.upr.date <- as.Date(newdat$pred.upr, origin = "2010-01-01")
df.mod.all.windows.start$b.start.date <- as.Date(df.mod.all.windows.start$b.start, origin = "2010-01-01")


Sys.setlocale("LC_TIME", "English")

# create plot with relative axis (by substracting the mean)

newdat$pred.relative <- as.numeric(newdat$pred.median.date - mean(newdat$pred.median.date))
newdat$pred.lwr.relative <- as.numeric(newdat$pred.lwr.date - mean(newdat$pred.median.date))
newdat$pred.upr.relative <- as.numeric(newdat$pred.upr.date - mean(newdat$pred.median.date))

df.mod.all.windows.start$b.start.date.relative <-  as.numeric(df.mod.all.windows.start$b.start.date - mean(df.mod.all.windows.start$b.start.date))


newdat$temp.relative <- newdat$temp.mean - mean(df.mod.all.windows.start$temp.mean)
df.mod.all.windows.start$temp.relative <- df.mod.all.windows.start$temp.mean - mean(df.mod.all.windows.start$temp.mean)

vec.breaks <- seq(-45, 60, by = 15)


ggplot.temp.prec.relative <- ggplot(newdat, aes(x=temp.relative, y = pred.relative, col = as.factor(prec.doy)))+
  geom_hline(aes(yintercept = 0), linetype = "dashed", col = "grey20")+
  geom_point(data = df.mod.all.windows.start, aes(temp.relative, b.start.date.relative), col = "darkgrey")+
  geom_line()+
  geom_ribbon(aes(ymin = pred.lwr.relative, ymax = newdat$pred.upr.relative, fill = as.factor(prec.doy)), alpha = 0.3, col = NA)+
  ylab("change in start of the breeding season (days)") + xlab("deviation from mean temperature °C \n (May 4 to June 12)")+
  xlim(-4.5,3)+
  scale_y_continuous(breaks = vec.breaks, limits = c(-45,68))+
  theme_classic()+ theme(legend.position = "top") + 
  labs(col = "number of precipitation days \n (April 18 - April 28)", 
       fill = "number of precipitation days \n (April 18 - April 28)")+ 
  theme(axis.text = element_text(size = 14), 
        axis.title =  element_text(size = 14),
        legend.title=element_text(size=14),
        legend.text = element_text(size=12))


pdf("temp.prec.relative.pdf")
print(ggplot.temp.prec.relative)     
dev.off() 

warning("Hey there! Thank you for having such a close look at the code :-)")

### END OF THE BREEDING SEASON ### Epred.medianND OF THE BREEDING SEASON 
### create a dataframe with evironemental variables for the time window with the max positive correlation coefficient

# time window for temperature

start.temp <- df.combinations.end$start[df.combinations.end$start.pcor.temp ==
                                          min(df.combinations.end$start.pcor.temp)]
end.temp <- df.combinations.end$end[df.combinations.end$start.pcor.temp == 
                                      min(df.combinations.end$start.pcor.temp)]


vec.doy.temp <- paste("doy", c(start.temp, end.temp) , sep = "" )


ind.start.temp <- which(colnames(df.temp.year) == vec.doy.temp[1])
ind.end.temp <- which(colnames(df.temp.year) == vec.doy.temp[2])

mean.temp <- rowMeans(df.temp.year[ind.start.temp:ind.end.temp])

# time window for precipitation

start.prec <- df.combinations.end$start[df.combinations.end$start.pcor.prec ==
                                          min(df.combinations.end$start.pcor.prec)]
end.prec <- df.combinations.end$end[df.combinations.end$start.pcor.prec ==
                                      min(df.combinations.end$start.pcor.prec)]


vec.doy.prec <- paste("doy", c(start.prec, end.prec) , sep = "" )


ind.start.prec<- which(colnames(df.prec.year) == vec.doy.prec[1])
ind.end.prec <- which(colnames(df.prec.year) == vec.doy.prec[2])

prec.doy <- rowSums(df.prec.year[ind.start.prec:ind.end.prec]>0, na.rm = T)
snow.day <- rowSums((df.prec.year[ind.start.prec:ind.end.prec]>0) & (df.temp.year[ind.start.prec:ind.end.prec]<0), na.rm = T)

  
df.mod.all.windows.end <- data.frame(b.length = df.SF$length.mean,
                     b.start = df.SF$start,
                     b.end = df.SF$end,
                       b.intensity = df.SF$intensity,
                       prec.doy = prec.doy,
                     snow.day = snow.day,
                       temp.mean = mean.temp,
                     snow.timing = df.SF$snow50,
                       melt.duration = df.SF$melt.duration,
                       year = df.SF$year,
                       breg = df.SF$bioreg)


mod.end <- brm(b.end ~ temp.mean  + snow.timing + melt.duration  + (1|year) + (1|breg), data = df.mod.all.windows.end)
summary(mod.end)
conditional_effects(mod.end)

#saveRDS(mod.end, file = "mod.end.rds")
mod.end <- readRDS("mod.end.rds")

newdat <- expand.grid(melt.duration = mean(df.mod.all.windows.end$melt.duration),
                     snow.timing = mean(df.mod.all.windows.end$snow.timing),
                     temp.mean = seq(min(df.mod.all.windows.end$temp.mean),
                                     to = max(df.mod.all.windows.end$temp.mean),
                                     length.out =100))# draw samples from the posterior distribution based on the newly created data
pred <- posterior_epred(mod.end, newdata = newdat, re_formula = NA)
# calculate the median and uncertainty intervall from the posterior draws
newdat$pred.median <- apply(pred, 2, mean)
newdat$pred.lwr <- apply(pred, 2, quantile, 0.025)
newdat$pred.upr <- apply(pred, 2, quantile, 0.975)



newdat$pred.median.date <- as.Date(newdat$pred.median, origin = "2010-01-01")
newdat$pred.lwr.date <- as.Date(newdat$pred.lwr, origin = "2010-01-01")
newdat$pred.upr.date <- as.Date(newdat$pred.upr, origin = "2010-01-01")
df.mod.all.windows.end$b.end.date <- as.Date(df.mod.all.windows.end$b.end, origin = "2010-01-01")


# create plot with relative axis (by substracting the mean)

newdat$pred.relative <- newdat$pred.median.date - mean(newdat$pred.median.date)
newdat$pred.lwr.relative <- newdat$pred.lwr.date - mean(newdat$pred.median.date)
newdat$pred.upr.relative <- newdat$pred.upr.date - mean(newdat$pred.median.date)

df.mod.all.windows.end$b.end.date.relative <-  df.mod.all.windows.end$b.end.date - mean(df.mod.all.windows.end$b.end.date)


newdat$temp.relative <- newdat$temp.mean - mean(df.mod.all.windows.end$temp.mean)
df.mod.all.windows.end$temp.relative <- df.mod.all.windows.end$temp.mean - mean(df.mod.all.windows.end$temp.mean)

vec.breaks <- seq(-45, 65, by = 15)


ggplot.temp.prec.relative <- ggplot(newdat, aes(x=temp.relative, y = pred.relative))+
  geom_hline(aes(yintercept = 0), linetype = "dashed", col = "grey20")+
  geom_point(data = df.mod.all.windows.end, aes(temp.relative, b.end.date.relative), col = "darkgrey")+
  geom_line(col = "#607298")+
  geom_ribbon(aes(ymin = pred.lwr.relative, ymax = newdat$pred.upr.relative), alpha = 0.45, col = NA, fill = "#607298")+
  ylab("change in end of the breeding season (days)") + xlab("deviation from mean temperature °C \n (June 14 to June 30)")+
  xlim(-4.5,4.2)+ scale_y_continuous(breaks = vec.breaks,  limits = c(-42,65))+
  theme_classic()+ theme(legend.position = "top") +
  theme(axis.text = element_text(size = 14), 
        axis.title =  element_text(size = 14))


pdf("temp.prec.relative.end.pdf")
print(ggplot.temp.prec.relative)     
dev.off() 


### create a dataframe with evironemental variables for the time window with the max negative (e.g min) correlation coefficient

# merge data sets with critical time windows for start and end of the breeding season

df.length <- merge(df.mod.all.windows.start, df.mod.all.windows.end, 
              by = c("year", "breg"))

df.length <- select(df.length, c(-snow.timing.y,
                            -melt.duration.y,
                            -b.length.y,
                            -b.start.y,
                            -b.end.y))



mod.duration <- brm(b.length.x ~ temp.mean.x * prec.doy.x + temp.mean.y + (1|year) + (1|breg), data = df.length)
summary(mod.duration)

#saveRDSmod.duration#saveRDS(mod.duration, file = "mod_duration.rds")


newdat <- expand.grid(temp.mean.x = seq(min(df.mod.all.windows.start$temp.mean), 
                                        max(df.mod.all.windows.start$temp.mean),
                                        length.out = 100),
                      prec.doy.x = mean(df.mod.all.windows.start$prec.doy),
                      temp.mean.y = mean(df.mod.all.windows.end$temp.mean),
                      snow.timing.x = mean(df.mod.all.windows.end$snow.timing),
                      melt.duration.x = mean(df.mod.all.windows.end$melt.duration))
# draw samples from the posterior distribution based on the newly created data
pred <- posterior_epred(mod.duration, newdata = newdat, re_formula = NA)


# calculate the median and uncertainty intervall from the posterior draws
newdat$pred.median <- apply(pred, 2, mean)
newdat$pred.lwr <- apply(pred, 2, quantile, 0.025)
newdat$pred.upr <- apply(pred, 2, quantile, 0.975)

newdat$pred.relative <- newdat$pred.median - mean(newdat$pred.median)
newdat$pred.lwr.relative <- newdat$pred.lwr - mean(newdat$pred.median)
newdat$pred.upr.relative <- newdat$pred.upr - mean(newdat$pred.median)



newdat2 <- expand.grid(temp.mean.x = mean(df.mod.all.windows.start$temp.mean),
                      prec.doy.x = mean(df.mod.all.windows.start$prec.doy),
                      temp.mean.y = seq(min(df.mod.all.windows.end$temp.mean), 
                                        max(df.mod.all.windows.end$temp.mean),
                                        length.out = 100),
                      snow.timing.x = mean(df.mod.all.windows.end$snow.timing),
                      melt.duration.x = mean(df.mod.all.windows.end$melt.duration))

# draw samples from the posterior distribution based on the newly created data
pred <- posterior_epred(mod.duration, newdata = newdat2, re_formula = NA)


# calculate the median and uncertainty intervall from the posterior draws
newdat2$pred.median <- apply(pred, 2, mean)
newdat2$pred.lwr <- apply(pred, 2, quantile, 0.025)
newdat2$pred.upr <- apply(pred, 2, quantile, 0.975)

newdat2$pred.relative <- newdat2$pred.median - mean(newdat2$pred.median)
newdat2$pred.lwr.relative <- newdat2$pred.lwr - mean(newdat2$pred.median)
newdat2$pred.upr.relative <- newdat2$pred.upr - mean(newdat2$pred.median)
newdat2$temp.y.relative <- newdat2$temp.mean.y - mean(df.length$temp.mean.y)
newdat2$temp.x.relative <- newdat2$temp.y.relative



df.length$b.length.relative <-  df.length$b.length.x - mean(df.length$b.length.x)
df.length$temp.mean.x.relative <- df.length$temp.mean.x - mean(df.length$temp.mean.x)
df.length$temp.mean.y.relative <- df.length$temp.mean.y - mean(df.length$temp.mean.y)



newdat$temp.x.relative <- newdat$temp.mean.x - mean(df.length$temp.mean.x)

ggplot.duration.relative <- ggplot(data = newdat, aes(x=temp.x.relative,
                          y=pred.relative))+
  geom_point(data = df.length, aes(x = temp.mean.x.relative, 
                                   y = b.length.relative), col = "#E2D1F9", alpha = 0.8)+
  geom_point(data = df.length, aes(x = temp.mean.y.relative, 
                                   y = b.length.relative), col = "#317773", alpha = 0.5)+
  geom_hline(aes(yintercept = 0), linetype = "dashed", col = "grey20")+
  
  geom_line(col = "#E2D1F9")+
  geom_ribbon(aes(ymin = pred.lwr.relative, ymax = newdat$pred.upr.relative), alpha = 0.5, col = NA, fill = "#E2D1F9")+
  
  geom_line(data = newdat2, aes(x= temp.y.relative,
                                y = pred.relative), col = "#317773", alpha = 0.75)+
  geom_ribbon(data = newdat2, aes(ymin = pred.lwr.relative, ymax = pred.upr.relative), alpha = 0.45, col = NA, fill = "#317773")+
  
  ylab("change in duration of breeding season (days)") + xlab("deviation from mean temperature °C")+
  xlim(-5,4.5)+
  theme_classic()+ theme(legend.position = "top") +
  theme(axis.text = element_text(size = 14), 
        axis.title =  element_text(size = 14))

pdf("duration.relative.pdf")
print(ggplot.duration.relative)     
dev.off() 


```

