# 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 = "." )

brainsize <- read_excel("brainSize.xlsx") 
  # data frame with information on body size and brain size from five sources
names.birds <- brainsize$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 brain size data
scr <- match(names.recorded, names.birds)
scr <- !is.na(scr)
names.recorded <- names.recorded[scr]
data.birdsC <- data.birds # object for saving species with complete data
data.birdsC <- data.birdsC[scr,]
data.birdsC <- data.birdsC[,-c(1,2)]

scr <- match(names.recorded, names.birds)
brainsize<- brainsize[scr,]
names.birds <- names.birds[scr]
  # species list 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 brain size for assemblages

# Objects ending with W are objects for brain sizes
# weighted by individuals

# meanResid = mean brain size corrected for body
#             size using allometric equation

richness <- apply(data.birds[,3:750] > 0, 2, sum)
meanResid <- numeric()
meanResidW <- numeric()
richnessWithData <- numeric()

for (i in (1: ncol(data.birdsC))) {
  scrD <- data.birdsC[,i] > 0
  scr <- ((data.birdsC[,i] > 0) * brainsize$residBrain)[which(scrD == T)]
  
  richnessWithData[i] <- length(scr) 
  meanResid[i] <- mean(scr)
  
  meanResidW[i] <- sum((data.birdsC[,i] * 
                   brainsize$residBrain)[which(scrD == T)])/
                   sum(data.birdsC[,i])
}
percWithData <- richnessWithData/richness

## plotting histogram of data completeness for individual assemblages
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)

hist(percWithData, main = "",
     xlab = "Completeness [% of richness]", 
     ylab = "Number of data sets",
     axes = F)
axis(1, at = c(0.8, 0.85, 0.9, 0.95, 1),
     labels = c( 80, 85, 90, 95, 100), lwd = 2)
axis(2, lwd =2)
box(lwd =2)
par(par.sav)

min(percWithData)
sum(percWithData > 0.9)
sum(percWithData > 0.9)/ncol(data.birdsC)


#########################################################
# 3a Plot of residuals from the regression 
#    log10 (brainsize) ~ log10(meanMass) using all data
#    available for the species ecoreded in Upper Franconia
#    against days starting January 01, 2010

#    red line: rolling mean using 3 brain size data

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, meanResid, axes = F,
     xlim = c(0, 5850),
     xlab = "Date - starting January 01, 2010",
     ylab = "Mean brain size - species",
     pch = 21, col = "darkgrey", bg = "darkgrey")
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, lwd =2)
box(lwd =2)
abline(v=(begin.years - dates("12/31/2009")), col = "grey")
rollLine <- rollmean(cbind(daySince2010, meanResid), k= 3)
lines(rollLine[,1], rollLine[,2], col = "red", lwd = 2)
par(par.sav) 

######################################################################
# 3b Plot of mean brain size weighted by the number of recorded
#    individuals against 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, meanResidW, axes = F,
     xlim = c(0, 5850),
     xlab = "Date - starting January 01, 2010",
     ylab = "Mean brain size - individuals",
     pch = 21, col = "darkgrey", bg = "darkgrey")
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, lwd =2)
box(lwd = 2)
abline(v=(begin.years - dates("12/31/2009")), col = "grey")

rollLine <- rollmean(cbind(daySince2010, meanResid), k= 3)
lines(rollLine[,1], rollLine[,2], col = "red", lwd = 2)
par(par.sav)

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

mean.D <- data.frame(date = daySince2010,
                     brain = meanResid)
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),
    #    mfcol = c(2,1),
    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.7),
     axes = F)
ticksx<-c(0, 182, 364, 500, 1000)
ticksy<-c(0, 0.2, 0.4, 0.6)
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)