############################################
# 1 Preparing data
# cleaning the environment
# loading packages

rm(list=ls(all=TRUE))
par.sav <- par(no.readonly = T)

library(readxl)    # reading excel files
library(chron)     # for handling dates
library(lubridate) # for calculating day of the year
library(lomb)      # for fitting a Lomb-Scargle periodogram
library(zoo)       # for calculating rolling mean

# reading data, matching names of the two data files
data <- read_excel("transectCountsDeps2010_2025.xlsx") 
data <-as.data.frame(data)
  ## basic data set of bird counts along one transect
data.birds <- data[,-1] 
  ## excluding the first column
names.recorded <- paste(data.birds$genus, data.birds$species, sep = "." )

bodyMass <- read_excel("bodyMass.xlsx")
  ## data frame with information on body size 
names.birds <- bodyMass$species

# excluding species with few records: n.records > n.limit
n.records <- rowSums(data.birds[,3:ncol(data.birds)])
  ## set limit for excluding birds with records of
  ## <= individuals recorded across all counts
n.limit <- 0
  ## at present no species excluded
keep.species <- n.records > n.limit
data.birds <- data.birds[keep.species,]
names.recorded <- names.recorded[keep.species]

# excluding species with no body size data
scr <- match(names.recorded, names.birds)
scr <- !is.na(scr)
names.recorded <- names.recorded[scr]
data.birdsC <- data.birds  # creating object for saving species with body mass data
data.birdsC <- data.birdsC[scr,]
data.birdsC <- data.birdsC[,-c(1,2)]

scr <- match(names.recorded, names.birds)
bodyMass<- bodyMass[scr,]
names.birds <- names.birds[scr]
  ## species in the two data files have the same ordering

# preparing tick information for the plots
begin.years <- dates(c("01/01/2010", "01/01/2011",
                       "01/01/2012", "01/01/2013",
                       "01/01/2014", "01/01/2015",
                       "01/01/2016", "01/01/2017",
                       "01/01/2018", "01/01/2019",
                       "01/01/2020", "01/01/2021",
                       "01/01/2022", "01/01/2023",
                       "01/01/2024", "01/01/2025",
                       "01/01/2026"))  

begin.years.long <- dates(c("01/01/10",
                            "01/01/14",
                            "01/01/18",
                            "01/01/22",
                            "01/01/26"))

begin.years.short <- setdiff(begin.years, begin.years.long)
  ## "short" and "long" refer to the tick length

year <- year(as.Date(chron(as.character(
  names(data.birdsC)),format=c(dates="d.m.y"))))
dayY <- yday(as.Date(chron(as.character(
  names(data.birdsC)),format=c(dates="d.m.y"))))
  ## days within a year
daySince2010 <- as.numeric(chron(as.character(names(data.birdsC)),
                                 format=c(dates="d.m.y"))
                           - dates("12/31/2009"))
  ## days since beginning of 2010

###################################################
# 2 Calculation of mean body size for assemblages

# Objects ending with W are objects for body sizes
# weighted by individuals
# for simplicity called "calculated across individuals"

# body size log10-transformed

meanbody <- numeric()
meanbodyW <- numeric()
for (i in (1: ncol(data.birdsC))) {
  scrD <- data.birdsC[,i] > 0
  scr <- ((data.birdsC[,i] > 0) * log10(bodyMass$meanMass))[which(scrD == T)]
  meanbody[i] <- mean(scr, na.rm =T)
  
  meanbodyW[i] <- sum((data.birdsC[,i] * 
                       log10(bodyMass$meanMass))[which(scrD == T)], na.rm = T)/
                       sum(data.birdsC[,i], na.rm = T)
}
##########################################################
# 3a Plot of community mean of body size across species
#    versus date starting January 01, 2010

par(mar = c(5, 6, 1, 1), mai = c(1, 1, 0.2, 0.2),
    lwd = 2, las = 1, cex = 0.9, 
    cex.axis = 1, cex.lab = 1.5)

plot(daySince2010, meanbody, axes = F,
     xlim = c(0, 5850),
     xlab = "Date - start January 01, 2010",
     ylab = "Mean body mass - species",
     pch = 21, col = "darkgrey", bg = "darkgrey")
yticks <- c(50, 100, 200, 400)
axis(1, at = (begin.years.long - dates("12/31/2009")), 
     labels = begin.years.long, lwd = 2, tck = -0.07)
axis(1, at = (begin.years.short - dates("12/31/2009")), 
     labels = rep("", length(begin.years.short)), lwd = 2,
     tck = -0.04)
axis(2, at = log10(yticks), labels = yticks, lwd =2)
box(lwd =2)
abline(v=(begin.years - dates("12/31/2009")), col = "grey")
rollLine <- rollmean(cbind(daySince2010, meanbody), k= 3)
lines(rollLine[,1], rollLine[,2], col = "red", lwd = 2)
par(par.sav)  

##########################################################
# 3b Plot of community mean of body mass weighted by
#    the number of recorded individuals

par(mar = c(5, 6, 1, 1), mai = c(1, 1, 0.2, 0.2),
    lwd = 2, las = 1, cex = 0.9, 
    cex.axis = 1, cex.lab = 1.5)

plot(daySince2010, meanbodyW, axes = F,
     xlim = c(0, 5850),
     xlab = "Date - start January 1, 2010",
     ylab = "Mean body mass - individuals",
     pch = 21, col = "darkgrey", bg = "darkgrey")
yticks <- c(50, 100, 200, 400)
axis(1, at = (begin.years.long - dates("12/31/2009")), 
     labels = begin.years.long, lwd = 2, tck = -0.07)
axis(1, at = (begin.years.short - dates("12/31/2009")), 
     labels = rep("", length(begin.years.short)), lwd = 2,
     tck = -0.04)
axis(2, at = log10(yticks), labels = yticks, lwd =2)
box(lwd =2)
abline(v=(begin.years - dates("12/31/2009")), col = "grey")
rollLine <- rollmean(cbind(daySince2010, meanbodyW), k= 3)
lines(rollLine[,1], rollLine[,2], col = "red", lwd = 2)
par(par.sav)  

#######################################################
# 3c Lomb-Scargle periodogram for mean size across
#    species using package lomb 

mean.D <- data.frame(date = daySince2010,
                     size = meanbody)
res.lsp <- lsp(mean.D, 
               from = 20, to = 1000,
               ofac = 5, type = "period", plot = F)
summary(res.lsp)  
(peaksData <- getpeaks(res.lsp, 2, F))

par(mar = c(5, 6, 1, 1), mai = c(1, 1, 0.2, 0.2),
    lwd = 2, las = 1,cex = 0.9, 
    cex.axis = 1, cex.lab = 1.5)

plot(res.lsp$scanned, 
     res.lsp$power,
     type = "l",
     xlab="period [days]", ylab="Normalized power - species",
     xlim = c(0, 1000), ylim = c(0, 0.3),
     axes = F)
ticksx<-c(0, 183.4, 364.5, 500, 1000)
ticksy<-c(0, 0.1, 0.2, 0.3)
axis(1, at=ticksx, labels=ticksx, lwd.ticks=2)
axis(2, at=ticksy, labels=ticksy, lwd.ticks=2)
abline(h = res.lsp$sig.level, col = "red")
text(600,res.lsp$sig.level + 0.08,
     labels = " p < 0.01",
     col = "red", cex = 1)
box()
par(par.sav)