---
title: "GAMMS for breeding phenology based on citizen science data"
author: "Carole Niffenegger"
date: "`r Sys.Date()`"
output: html_document
---

Current R version: 4.2.2

## Load packages

```{r}

library(ggplot2)
library(dplyr)
library(birdring)
library(parallel)
library(readr)
library(rstanarm)
library(gridExtra)
library(mgcv)
library(brms)
library(bayesplot)
library(ggpubr)
library(GGally)


```

## Data analysis
### GAM with day of year

[Reference]{.underline}: Pedersen EJ, Miller DL, Simpson GL, Ross N. 2019. Hierarchical generalized additive models in ecology: an introduction with mgcv. PeerJ 7:e6876 <https://doi.org/10.7717/peerj.6876>

```{r, eval = FALSE}

SF_breeding <- read_delim("SF_breeding.csv",  delim = ";", escape_double = FALSE, trim_ws = TRUE)

# add a column "day of year"

SF_breeding$doy <- dayofyear(SF_breeding$monat, SF_breeding$tag, SF_breeding$jahr)

# this model fits one GAM per year (all regions combined)

# The following loop calculates a simple GAM with the binary variable "BreedingInd" as response variable and 
# day of year (e.g. date) as explanatory variable
# I fit one GAM per year using the brm function from the brms package

# I use s() for the smoother with cyclic cubic regression splines (specified as bs = "cc"). The start and end of cc smoother 
# are constrained to match in values and first derivatives. This seems usefull in this case as I want the breeding intensity
# to start and end at the same level (ideally at zero). 

warning("Adjust number of cores")

cores <- detectCores()-2

# format "Jahr" as factor
SF_breeding$jahr <- as.factor(SF_breeding$jahr)

# Estimated run time: 80 minutes

for(for.year in 1999:2022){
  SF <- SF_breeding[SF_breeding$jahr == for.year,]
  file.name <- paste("mod", for.year, "univariate.rda", sep = "")
  mod <- brm(breedingInd ~ s(doy, bs = "cc", k = 5) , family = "bernoulli", data = SF,
             cores = cores, iter = 40000, thin = 4, control = list(adapt_delta = 0.9), chains = 4)
  saveRDS(mod, file = file.name)
}


```

###### MCMC diagnostics

```{r, eval = FALSE}

# I did not find a good documentation if or how to do MCMC diagnostics for Generalized additive models. Therefore I used the standard plot as I would for a linear regression.

# since I fitted one GAM per year, I have many diagnostic plots
# The following loop produces traceplots, acf plots and R hat values for every model separately

# create list of all models
mod.list <- list.files(pattern = "univariate.rda")

rhat.univariate <- as.data.frame(matrix(nrow= 0, ncol = 8))
column.names <- c("model", "Rhat_intercept", "Rhat_sds_sdoy_1", "Rhat_s_sdoy_1[1]","Rhat_s_sdoy_1[2]",
                  "Rhat_s_sdoy_1[3]", "Rhat_lprior", "Rhat_lp__")
colnames(rhat.univariate) <- column.names


for (for.diagnostics in mod.list) {
  mod <- readRDS(for.diagnostics)
  
  # produce trace plot and save as pdf
  pdf(paste("traceplot_", for.diagnostics, ".pdf", sep = ""))
  print(mcmc_trace(mod) + ggtitle(for.diagnostics)+ theme(plot.title = element_text(hjust = 0.5)))
  dev.off()
  
  # produce acf plot and save as pdf
  
  pdf(paste("acf_", for.diagnostics, ".pdf", sep = ""))
  print(mcmc_acf_bar(mod)+ ggtitle(for.diagnostics)+ theme(plot.title = element_text(hjust = 0.5)))
  dev.off()
  
  rhat.univariate[(nrow(rhat.univariate)+1),] <- c(for.diagnostics, rhat(mod))
  
}

sum(rhat.univariate[,2:8] > 1.01) # 0

# Overall, the mcmc diagnostics look ok


```

###### Visualize results

```{r}

mod.list <- list.files(pattern = "univariate.rda")
plot_list <- list()

for (i in 1:length(mod.list)) {
  for.plot <- mod.list[i]
  # read model
  mod <- readRDS(for.plot)
  # produce plot wiht the effect of date on the breeding index based on the 
  # posterior distribution
  effect_plot <- conditional_effects(mod, prob = 0.95, robust = FALSE)
  # prob : probability covered in CrI (0.95 is the default value)
  # robust = FALSE: the regression lines is the mean rather than the median
  # extract the year for plot title
  title <- substr(for.plot,4,7)
  # modify plot and save in a list
  plot_list[[i]] <- plot(effect_plot, plot = FALSE)[[1]] +
  ylim(0,1)+
  theme_bw()+
  ggtitle(title)
  
}

grid.arrange(grobs=plot_list,ncol=6)

# The model in 2000 seems problematic. Base level is too high.
# Otherwise, the results seems plausible

```

##### Calculate length of the breeding season

```{r}

# Here, I calculate the length of the breeding season for every year
# based on the previously fitted models
# By drawing samples from the posterior distribution, I can additionally # calculate credible intervals for the lenght of the breeding season

# NOTE: I need to somehow define, what I consider the lenght of the breeding season. 
# There are potentially many options - some more elegant than others ;-)
# For simplicity I use the following definition: breeding season starts when the breeding
# index exceeds 0.05 for the first time and ends when it is below 0.05 again.
# Again, it would be worth thinking about a more sophisticated / less random definition,
# for example based on the credible interval. But for now, this should give some first insights

# create a function that calculates the lenght of the breeding season based on the definition above
# secondly, the function returns the first day of the breeding season (start)
breeding_length <- function(x){
  start <- min(which(x > 0.05))
  end <- max(which(x > 0.05))
  breeding_length <- end-start
  return(c(breeding_length, start,end))
}


df.length.breeding <- data.frame(year = as.numeric(),
                                 mean = as.numeric(),
                                 lwr = as.numeric(),
                                 upr = as.numeric(),
                                 start = as.numeric(),
                                 end = as.numeric())


mod.list <- list.files( pattern = "univariate.rda")

for (for.dat in mod.list) {
  # import model
  mod <-readRDS(for.dat)
  # predict from the posterior predictive distribution
  # by using the posterior_epred() function, I only incorporate uncertainty in the 
  # expected value but ignore the uncertainty of the sampling process (e.g. the residuals)
  dat <- posterior_epred(mod, newdata = data.frame(doy = 91:273)) # doy from 1st of   April until 30th of September
  # calculate the lenght of the breeding seson for every draw using the function defined before
  fit.lenght <- apply(dat, 1, breeding_length)
  year <- as.numeric(substr(for.dat,4,7))
  # mean lenght
  mean <- mean(fit.lenght[1,])
  # 95% uncertaintiy interval
  quantiles <- as.numeric(quantile(fit.lenght[1,], probs = c(0.025,0.975)))
  start_day <- mean(fit.lenght[2,])+91 # here I add 91 because day 1 in the prediction corresponds to day 91
  end_day <- mean(fit.lenght[3,])+91
  df.length.breeding[(nrow(df.length.breeding)+1),] <- c(year, mean, quantiles, start_day, end_day)
}

write.csv(df.length.breeding, "SF_length.csv")

# summary statistics 
mean(df.length.breeding$mean)
range(df.length.breeding$mean)

# plot the lenght of the breeding season per year

# For the moment, I remove the year 2000 as there were some problems with the model
# the baseline was too high never reaching 0 even during times of the year when there was certainly no breeding

df.length.breeding[2, 2:5] <- NA

ggplot(df.length.breeding, aes(x = year, y = mean, ymin = lwr, ymax = upr))+ 
  geom_pointrange()+
  ylim(0,250)+
  ylab("Length (days)")+
  theme_bw()+
  geom_hline(yintercept = mean(df.length.breeding$mean, na.rm = T),
             col = "red", linetype = "dashed")+
  geom_text(x=2017, y=242, label="range (mean) = 70 - 135 days")


```


### GAM with biogeographic region

##### Fit the model

```{r}

SF_breeding$DERegionNa <- as.factor(SF_breeding$DERegionNa)

cores <- detectCores() - 2

for (for.year in 1999:2022) {
  # subset data (one year)
  for.SF <- SF_breeding[SF_breeding$jahr == for.year,]
  # fit GAM
  for.mod <- brm(breedingInd ~ s(doy, bs = "cc", k = 5) + 
                   s(doy, by  = DERegionNa, bs = "cc", k = 5, m = 1) + 
                   s(DERegionNa, bs = "re") , chains = 4, control = list(adapt_delta = 0.90),
                 iter = 80000, thin = 4,
           data = for.SF , family = "bernoulli", cores = cores)
  # create file name to save model
  for.file.name <- paste("mod_Breg", for.year, ".rda", sep = "")
  saveRDS(for.mod, file = for.file.name)
}

# based on the output, I will only use years 2006-2021 for further analysis
# the number of observations for this time window is


sum(SF_breeding$jahr >= 2006 & SF_breeding$jahr <= 2021)
SF_breeding0621 <- SF_breeding[SF_breeding$jahr >= 2006 & SF_breeding$jahr <= 2021,]

table(SF_breeding0621$breedingInd)

```


##### Calculate length of breeding season

```{r}

breeding_length <- function(x){
  start <- min(which(x > 0.05))
  end <- max(which(x > 0.05))
  breeding_length <- end-start
  return(c(breeding_length, start,end))
}



df.length.breeding.Breg <- data.frame(year = as.numeric(),
                                 Breg = as.character(),
                                 mean = as.numeric(),
                                 lwr = as.numeric(),
                                 upr = as.numeric(),
                                 start = as.numeric(),
                                 end = as.numeric())

mod.list4 <- list.files(pattern = "mod_Breg")



for (for.dat in mod.list4) {
  # import model
  mod <- readRDS(for.dat)
  # predict from the posterior predictive distribution
  # by using the posterior_epred() function, I only incorporate uncertainty in the 
  # expected value but ignore the uncertainty of the sampling process (e.g. the residuals)
  newdat <- expand.grid(doy = 91:273, DERegionNa = levels(as.factor(SF_breeding$DERegionNa)))
  dat <- posterior_epred(mod, newdata = newdat) # doy from 1st of April until 30th of September
  # split data set according the elevation
  datO <- dat[,which(newdat$DERegionNa == "Östliche Zentralalpen")]
  datS <- dat[,which(newdat$DERegionNa == "Alpensüdflanke")]
  datN <- dat[,which(newdat$DERegionNa == "Alpennordflanke")]
  datW <- dat[,which(newdat$DERegionNa == "Westliche Zentralalpen")]
  
  
  # calculate the lenght of the breeding season for every draw using the function defined before
  fit.lenghtO <- apply(datO, 1, breeding_length)
  year <- as.numeric(substr(for.dat,9,12))
  # mean lenght
  mean <- mean(fit.lenghtO[1,], na.rm = TRUE)
  # 95% uncertaintiy interval
  quantiles <- as.numeric(quantile(fit.lenghtO[1,], probs = c(0.025,0.975), na.rm = TRUE))
  start_day <- mean(fit.lenghtO[2,], na.rm = TRUE)+91 # here I add 91 because day 1 in the prediction corresponds to day 91
  end_day <- mean(fit.lenghtO[3,], na.rm = TRUE)+91 # here I add 91 because day 1 in the prediction corresponds to day 91
  
  df.length.breeding.Breg[(nrow(df.length.breeding.Breg)+1),] <- c(year, "Östliche Zentralalpen", mean, quantiles, start_day, end_day)
  
  
  # calculate the lenght of the breeding seson for every draw using the function defined before
  fit.lenghtS <- apply(datS, 1, breeding_length)
  year <- as.numeric(substr(for.dat,9,12))
  # mean lenght
  fit.lenghtS[1,][!is.finite(fit.lenghtS[1,])] <- NA
  mean <- mean(fit.lenghtS[1,], na.rm = TRUE)
  # 95% uncertaintiy interval
  quantiles <- as.numeric(quantile(fit.lenghtS[1,], probs = c(0.025,0.975), na.rm = TRUE))
  fit.lenghtS[2,][!is.finite(fit.lenghtS[2,])] <- NA
  start_day <- mean(fit.lenghtS[2,], na.rm = TRUE)+91 # here I add 91 because day 1 in the prediction corresponds to day 91
  fit.lenghtS[3,][!is.finite(fit.lenghtS[3,])] <- NA
  end_day <- mean(fit.lenghtS[3,], na.rm = TRUE)+91 # here I add 91 because day 1 in the prediction corresponds to day 91
  
  df.length.breeding.Breg[(nrow(df.length.breeding.Breg)+1),] <- c(year, "Alpensüdflanke", mean, quantiles, start_day, end_day)
  
  
  # calculate the lenght of the breeding seson for every draw using the function defined before
  fit.lenghtN <- apply(datN, 1, breeding_length)
  year <- as.numeric(substr(for.dat,9,12))
  # mean lenght
  fit.lenghtN[1,][!is.finite(fit.lenghtN[1,])] <- NA
  mean <- mean(fit.lenghtN[1,], na.rm = TRUE)
  # 95% uncertaintiy interval
  quantiles <- as.numeric(quantile(fit.lenghtN[1,], probs = c(0.025,0.975), na.rm = TRUE))
  fit.lenghtN[2,][!is.finite(fit.lenghtN[2,])] <- NA
  start_day <- mean(fit.lenghtN[2,], na.rm = TRUE)+91 # here I add 91 because day 1 in the prediction corresponds to day 91
  fit.lenghtN[3,][!is.finite(fit.lenghtN[3,])] <- NA
  end_day <- mean(fit.lenghtN[3,], na.rm = TRUE)+91 # here I add 91 because day 1 in the prediction corresponds to day 91
  
  df.length.breeding.Breg[(nrow(df.length.breeding.Breg)+1),] <- c(year, "Alpennordflanke", mean, quantiles, start_day, end_day)
  
 # calculate the lenght of the breeding seson for every draw using the function defined before
  fit.lenghtW <- apply(datW, 1, breeding_length)
  year <- as.numeric(substr(for.dat,9,12))
  # mean lenght
  fit.lenghtW[1,][!is.finite(fit.lenghtW[1,])] <- NA
  mean <- mean(fit.lenghtW[1,], na.rm = TRUE)
  # 95% uncertaintiy interval
  quantiles <- as.numeric(quantile(fit.lenghtW[1,], probs = c(0.025,0.975), na.rm = TRUE))
  fit.lenghtW[2,][!is.finite(fit.lenghtW[2,])] <- NA
  start_day <- mean(fit.lenghtW[2,], na.rm = TRUE)+91 # here I add 91 because day 1 in the prediction corresponds to day 91
  fit.lenghtW[3,][!is.finite(fit.lenghtW[3,])] <- NA
  end_day <- mean(fit.lenghtW[3,], na.rm = TRUE)+91 # here I add 91 because day 1 in the prediction corresponds to day 91
  
  df.length.breeding.Breg[(nrow(df.length.breeding.Breg)+1),] <- c(year, "Westliche Zentralalpen", mean, quantiles, start_day, end_day)
    
  
}

#write.csv(df.length.breeding.Breg, "SF_length.Breg.csv")


```

##### Test for observer effects

```{r}

# in this section, I will re-compute the duration of the breeding season while removing observations from the data set at random
# If the results remain relatively similary even if I remove data, I assume that the number of observations was sufficient to estimate the phenology curvs

for.year <- 2020
for.percentage <- 10

percentage.reduction <- c(2,5,10,15,20,30,50)

df.length.breeding.Breg.sens <- data.frame(year = as.numeric(),
                                 Breg = as.character(),
                                 mean = as.numeric(),
                                 lwr = as.numeric(),
                                 upr = as.numeric(),
                                 start = as.numeric(),
                                 end = as.numeric(),
                                 percentage = as.numeric())


for (for.year in 2006:2021) {
  for.SF <- SF_breeding[SF_breeding$jahr == for.year,]
  for (for.percentage in percentage.reduction) {
    n.red <- nrow(for.SF)*(for.percentage/100)
    for.SF.reduced <- for.SF[-sample(1:nrow(for.SF), n.red),]
    
    for.mod <- brm(breedingInd ~ s(doy, bs = "cc", k = 5) + 
                   s(doy, by  = DERegionNa, bs = "cc", k = 5, m = 1) + 
                   s(DERegionNa, bs = "re") , chains = 4, 
                   control = list(adapt_delta = 0.90),
           data = for.SF.reduced , family = "bernoulli", cores = cores)
    
    conditional_effects(for.mod)
    
    # expected value but ignore the uncertainty of the sampling process (e.g. the residuals)
  newdat <- expand.grid(doy = 91:273, DERegionNa = levels(as.factor(SF_breeding$DERegionNa)))
  dat <- posterior_epred(for.mod, newdata = newdat) # doy from 1st of April until 30th of September
  # split data set according the elevation
  datO <- dat[,which(newdat$DERegionNa == "Östliche Zentralalpen")]
  datS <- dat[,which(newdat$DERegionNa == "Alpensüdflanke")]
  datN <- dat[,which(newdat$DERegionNa == "Alpennordflanke")]
  datW <- dat[,which(newdat$DERegionNa == "Westliche Zentralalpen")]
  
  
  # calculate the lenght of the breeding season for every draw using the function defined before
  fit.lenghtO <- apply(datO, 1, breeding_length)
  year <- as.numeric(substr(for.dat,9,12))
  # mean lenght
  mean <- mean(fit.lenghtO[1,], na.rm = TRUE)
  # 95% uncertaintiy interval
  quantiles <- as.numeric(quantile(fit.lenghtO[1,], probs = c(0.025,0.975), na.rm = TRUE))
  start_day <- mean(fit.lenghtO[2,], na.rm = TRUE)+91 # here I add 91 because day 1 in the prediction corresponds to day 91
  end_day <- mean(fit.lenghtO[3,], na.rm = TRUE)+91 # here I add 91 because day 1 in the prediction corresponds to day 91
  
  df.length.breeding.Breg.sens[(nrow(df.length.breeding.Breg.sens)+1),] <- c(for.year, "Östliche Zentralalpen", mean, quantiles, start_day, end_day, for.percentage)
  
  
  # calculate the lenght of the breeding seson for every draw using the function defined before
  fit.lenghtS <- apply(datS, 1, breeding_length)
  year <- as.numeric(substr(for.dat,9,12))
  # mean lenght
  fit.lenghtS[1,][!is.finite(fit.lenghtS[1,])] <- NA
  mean <- mean(fit.lenghtS[1,], na.rm = TRUE)
  # 95% uncertaintiy interval
  quantiles <- as.numeric(quantile(fit.lenghtS[1,], probs = c(0.025,0.975), na.rm = TRUE))
  fit.lenghtS[2,][!is.finite(fit.lenghtS[2,])] <- NA
  start_day <- mean(fit.lenghtS[2,], na.rm = TRUE)+91 # here I add 91 because day 1 in the prediction corresponds to day 91
  fit.lenghtS[3,][!is.finite(fit.lenghtS[3,])] <- NA
  end_day <- mean(fit.lenghtS[3,], na.rm = TRUE)+91 # here I add 91 because day 1 in the prediction corresponds to day 91
  
  df.length.breeding.Breg.sens[(nrow(df.length.breeding.Breg.sens)+1),] <- c(for.year, "Alpensüdflanke", mean, quantiles, start_day, end_day, for.percentage)
  
  
  # calculate the lenght of the breeding seson for every draw using the function defined before
  fit.lenghtN <- apply(datN, 1, breeding_length)
  year <- as.numeric(substr(for.dat,9,12))
  # mean lenght
  fit.lenghtN[1,][!is.finite(fit.lenghtN[1,])] <- NA
  mean <- mean(fit.lenghtN[1,], na.rm = TRUE)
  # 95% uncertaintiy interval
  quantiles <- as.numeric(quantile(fit.lenghtN[1,], probs = c(0.025,0.975), na.rm = TRUE))
  fit.lenghtN[2,][!is.finite(fit.lenghtN[2,])] <- NA
  start_day <- mean(fit.lenghtN[2,], na.rm = TRUE)+91 # here I add 91 because day 1 in the prediction corresponds to day 91
  fit.lenghtN[3,][!is.finite(fit.lenghtN[3,])] <- NA
  end_day <- mean(fit.lenghtN[3,], na.rm = TRUE)+91 # here I add 91 because day 1 in the prediction corresponds to day 91
  
  df.length.breeding.Breg.sens[(nrow(df.length.breeding.Breg.sens)+1),] <- c(for.year, "Alpennordflanke", mean, quantiles, start_day, end_day, for.percentage)
  
 # calculate the lenght of the breeding seson for every draw using the function defined before
  fit.lenghtW <- apply(datW, 1, breeding_length)
  year <- as.numeric(substr(for.dat,9,12))
  # mean lenght
  fit.lenghtW[1,][!is.finite(fit.lenghtW[1,])] <- NA
  mean <- mean(fit.lenghtW[1,], na.rm = TRUE)
  # 95% uncertaintiy interval
  quantiles <- as.numeric(quantile(fit.lenghtW[1,], probs = c(0.025,0.975), na.rm = TRUE))
  fit.lenghtW[2,][!is.finite(fit.lenghtW[2,])] <- NA
  start_day <- mean(fit.lenghtW[2,], na.rm = TRUE)+91 # here I add 91 because day 1 in the prediction corresponds to day 91
  fit.lenghtW[3,][!is.finite(fit.lenghtW[3,])] <- NA
  end_day <- mean(fit.lenghtW[3,], na.rm = TRUE)+91 # here I add 91 because day 1 in the prediction corresponds to day 91
  
  df.length.breeding.Breg.sens[(nrow(df.length.breeding.Breg.sens)+1),] <- c(for.year, "Westliche Zentralalpen", mean, quantiles, start_day, end_day, for.percentage)
    
  
    
    
  }
  
}

write.csv(df.length.breeding.Breg.sens, "df.length.sensitiviy.csv")

```

##### Uncertainty of smoother

```{r}

library(brms)

mod.list5 <- list.files(pattern = "mod_Breg")
df.phenology <- read.csv("df.length.sensitiviy.csv")

# subset to years that were used in the analysis (2006-2021)

mod.list5 <- mod.list5[8:23]
df.phenology <- df.phenology[df.phenology$year >= 2006,]
df.phenology <- df.phenology[df.phenology$percentage == 5,] # because file contains several thresholds (which was used to test if the threshold influences the )

phenology <- data.frame(x = integer(),
                        year = integer(),
                        Breg = character(),
                        mean = numeric(),
                        lwr = numeric(),
                        upr = numeric(),
                        start = numeric(),
                        end = numeric(),
                        percentage = integer(),
                        uncertainty = numeric())

# create a loop that extracts the uncertainty interval for every year and region

for.mod <- mod.list5[3]

for (for.mod in mod.list5) {
  # load model
  
  mod <- readRDS(for.mod)
  for.year <- as.numeric(substr(for.mod,9,12))
 
  newdat <- expand.grid(doy = 91:273, DERegionNa = levels(as.factor(SF_breeding$DERegionNa)))
  dat <- posterior_epred(mod, newdata = newdat) 
  
  # split data set according the elevation
  datO <- dat[,which(newdat$DERegionNa == "Östliche Zentralalpen")]
  datS <- dat[,which(newdat$DERegionNa == "Alpensüdflanke")]
  datN <- dat[,which(newdat$DERegionNa == "Alpennordflanke")]
  datW <- dat[,which(newdat$DERegionNa == "Westliche Zentralalpen")]
   
  for.phenology <- df.phenology[df.phenology$year == for.year,]
  
  # extract sd of predictions for the START of breeding
  
  for.startO <- for.phenology$start[1]-91
  for.startS <- for.phenology$start[2]-91
  for.startN <- for.phenology$start[3]-91
  for.startW <- for.phenology$start[4]-91
  
  
sdO <- sd(datO[,for.startO])
sdS <- sd(datS[,for.startS])
sdN <- sd(datN[,for.startN])
sdE <- sd(datW[,for.startW])

  # extract sd of predictions for the END of breeding
  
  for.endO <- for.phenology$end[1]-91
  for.endS <- for.phenology$end[2]-91
  for.endN <- for.phenology$end[3]-91
  for.endW <- for.phenology$end[4]-91
  
  
sdO_end <- sd(datO[,for.endO])
sdS_end <- sd(datS[,for.endS])
sdN_end <- sd(datN[,for.endN])
sdE_end <- sd(datW[,for.endW])


for.phenology$start_uncertainty <- c(sdO, sdS, sdN, sdE)
for.phenology$end_uncertainty <- c(sdO_end, sdS_end, sdN_end, sdE_end)
  
phenology <- rbind(phenology, for.phenology)

  
}

write.csv(phenology, "uncertainty_start_end.csv")

plot(phenology$start, phenology$start_uncertainty)
cor(phenology$start, phenology$start_uncertainty)

plot(phenology$end, phenology$end_uncertainty)
cor(phenology$end, phenology$end_uncertainty)


# create plots about the correlation of uncertainty and the start and end date

pdf("uncertainty_start.pdf")
ggplot(data= phenology, aes(x=start, y = start_uncertainty))+
  geom_smooth(col = "red", size = 0.5)+  
  geom_point()+
  xlab("start day (day of year)")+ ylab("sd of smoother")+
  theme_classic()
dev.off()

pdf("uncertainty_end.pdf")
ggplot(data= phenology, aes(x=end, y = end_uncertainty))+
  geom_smooth(col = "red", size = 0.5)+ 
  geom_point()+
  xlab("end day (day of year)")+ ylab("sd of smoother")+
  theme_classic()
dev.off()

```

##### Semi-variograms

```{r}

library(gstat)
library(sp)

# produce semi-variograms of all models

mod.list <- list.files(pattern = ".rda")

# since I only analysed data starting in 2006, I remove all other models here

mod.list.short <- mod.list[8:24]


for.mod <- mod.list.short[2]

for (for.mod in mod.list.short) {
  
  # load model
  mod <- readRDS(for.mod)
  # extract year from model name
  for.year <- as.numeric(substr(for.mod, 9,12))
  # select SF data from year (to get coordinates of observations)
  for.SF <- SF_breeding[SF_breeding$jahr == for.year,]
  # create a data frame with model residuals and coordinates
  resdat <- data.frame(resid=resid(mod)[,1], x=for.SF$koordx,
                     y=for.SF$koordy)
  # convert coordinates in spatial object
  coordinates(resdat) = ~x+y
  # calculate semi-variances
  vario.mod <- variogram(resid~x+y, data = resdat, width = 2, cutoff = 20)
  # save plot
  file.name <- paste("variograms/variogram", for.year, "png",
                     sep = ".")
  png(filename=file.name)
  print(plot(vario.mod))
  dev.off()
  
  
}




```

##### ACF plot

```{r}

# this analysis can only be performed with the original data set containing the coordinates of each observation. To get access to the uncensored 
# data, please contact carole.niffenegger@vogelwarte.ch

# produce semi-variograms of all models

mod.list <- list.files(pattern = ".rda")

# since I only analysed data starting in 2006, I remove all other models here

mod.list.short <- mod.list[8:24]

for.mod <- mod.list.short[2]

for (for.mod in mod.list.short) {
  
  # load model
  mod <- readRDS(for.mod)
  # extract year from model name
  for.year <- as.numeric(substr(for.mod, 9,12))
  # save residuals
  residuals <- resid(mod)[,1]
  # save plot
  file.name <- paste("acf plots/ACF", for.year, "png",
                     sep = ".")
  png(filename=file.name)
  print(acf(residuals))
  dev.off()
  
  
}

```

