#--------------------------------------------------------------------------------------------------------------#
# function to estimate flowering overlap (with overlapEst) between all pairs of plant species within each year #
#--------------------------------------------------------------------------------------------------------------#

library(overlap)


## function that returns all overlap estimates (Dhat1, Dhat4, Dhat5) (see Ridout & Linkie 2009)

overlap<-function(mydata)
{
  results<-NULL
  mydata<-read.table(mydata,header=T,sep="\t", stringsAsFactors = FALSE)
  #select unique years
  y<-unique(mydata$Year)
  
  
  #for each year
  for(annee in y)
  {
    print("annee")
    print(annee)
    #select all lines in the table 'mydata' that correspond to the year annee
    subset_year<-mydata[which(mydata$Year == annee),]
    sp<-unique(subset_year$Species)#select species from the year "annee" (of the loop)
    espece_done<-NULL #vector that contains species already compared (to avoid replicates). It is empty, it will fill as comparisons occur
    #for every espece1
    for(espece in sp)
    {
      espece_done<-cbind(espece_done, espece)
      print("espece1")
      print(espece)
      #select all lines of the sub-table "annee" that correspond to the species specie
      #select the species to be compared with
      subset_specie<- subset_year[which(subset_year$Species == espece),] #subset of espece1 within a year
      #compare to all the others
      sp2<- sp[- which(sp %in% espece_done)] #sp2: all the other species that have not been compared yet
      if(length(sp2>1)) #if there is any left to compare
      {
        #for each species that has not been compared
        for(espece2 in sp[- which(sp %in% espece_done)])
        {
          print("espece2")
          print(espece2)
          print("ok")
          subset_specie2<- subset_year[which(subset_year$Species == espece2),] #species2 to be compared with specie1
          res<-as.matrix(overlapEst( subset_specie$doyRad,  subset_specie2$doyRad)) #result of the function overlapEst
          colnames(res)<-paste(paste(espece, espece2, sep="_and_"), annee, sep="-") #add column name with different separators
          results<-cbind(results, res) #add this result to the final table
        }
      }
    }
  }
  return(results)
}


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


## automatically select Dhat1 or Dhat4 based on species abundance (see Ridout & Linkie 2009)


overlap2<-function(mydata)
{
  results<-NULL
  mydata<-read.table(mydata,header=T,sep="\t", stringsAsFactors = FALSE)
  #select unique years
  y<-unique(mydata$Year)
  
  
  #for each year
  for(annee in y)
  {
    #print("annee")
    #print(annee)
    #select all lines in the table 'mydata' that correspond to the year annee
    subset_year<-mydata[which(mydata$Year == annee),]
    sp<-unique(subset_year$Species)#select species from the year "annee" (of the loop)
    espece_done<-NULL #vector that contains species already compared (to avoid replicates). It is empty, it will fill as comparisons occur
    #for every espece1
    for(espece in sp)
    {
      espece_done<-cbind(espece_done, espece)
      #print("espece1")
      #print(espece)
      #select all lines of the sub-table "annee" that correspond to the species specie
      #select the species to be compared with
      subset_specie<- subset_year[which(subset_year$Species == espece),] #subset of espece1 within a year
      #compare to all the others
      sp2<- sp[- which(sp %in% espece_done)] #sp2: all the other species that have not been compared yet
      if(length(sp2>1)) #if there is any left to compare
      {
        #for each species that has not been compared
        for(espece2 in sp[- which(sp %in% espece_done)])
        {
          #print("espece2")
          #print(espece2)
          #print("ok")
          subset_specie2<- subset_year[which(subset_year$Species == espece2),] #species2 to be compared with specie1

           if(dim(subset_specie)[1] > 50  & dim(subset_specie2)[1] > 50) 
             {
             #print("risotto")
             res1<-as.matrix(overlapEst( subset_specie$doyRad,  subset_specie2$doyRad, type = "Dhat4")) #result of the function overlapEst
           colnames(res1)<-paste(paste(espece, espece2, sep="_and_"), annee, sep="-") #add column name with different separators
           results<-cbind(results, res1) #add this result to the final table
           }
          
          else 
          {
           #print("pizza") 
            res2<-as.matrix(overlapEst( subset_specie$doyRad,  subset_specie2$doyRad, type = "Dhat1")) #result of the function overlapEst
          colnames(res2)<-paste(paste(espece, espece2, sep="_and_"), annee, sep="-") #add column name with different separators
          results<-cbind(results, res2) #add this result to the final tableresults<-cbind(results, res) #add this result to the final table
          }
          
          #results3<-cbind(results, results2) #add this result to the final tableresults<-cbind(results, res) #add this result to the final table
          
          
            }
      }
    }
  }
 rownames(results)<-"overlap_Est"
  #return(res1)
  return(results)
}



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

# -------------------------------------------------------------------------#
# Estimate flowering overlap within elevation band (i.e. within community) # 
# -------------------------------------------------------------------------#

## Load packages and dataset##
library(overlap)
library(tidyr)
library(tibble)
library(dplyr)
library(lubridate)



total <- read.table("total_rawdata.txt", sep = "\t", header = T)


#convert "Hike_date" to Julian day
total$date <- as.Date(total$Hike_date, format = "%d/%m/%Y") #Capital Y for 4-digit years
total$doy <- yday(total$date)  

total <- select(total, c("Mile", "Species", "Year", "doy"))

#-------------------------------------------------------------#
## convert doy to radians accounting for leap/non-leap years ##
#-------------------------------------------------------------#

#plants$leap <- (plants$doy / 366) * 360  #check for leap years transformation
#plants$nonleap <- (plants$doy / 365) * 360 #check for non-leap year transformation

total$is_leap <- as.character(leap_year(total$Year)) #check for leap/non-leap years

total$doyDeg <- ifelse(total$is_leap == TRUE, (total$doy / 366) * 360, (total$doy / 365) * 360) #convert doy in degrees accounting for leap years

total$doyRad <- total$doyDeg * pi / 180 #convert degrees in radians 

total<-select(total, c("Mile" ,"Species", "Year", "doyRad"))

#----------------------------------------------------------------------#
## exclude species in any year in which they appear less than 4 times ##
#----------------------------------------------------------------------#

total <- mutate(total,
                     spym = paste(Species, Year, Mile, sep = '-'))  #concatenate specie, year and mile in new column

final <- total %>% group_by(spym) %>% filter(n()>=4) #filter by new concatenated factors (species*year*mile) that appear >= 4 times 

sp_occurrence <- as.data.frame(table(sub$spym)) #frequence of every species in each year*mile where they occur
min(sp_occurrence$Freq)
max(sp_occurrence$Freq)

final<-final[,-5]

write.table(final, file = "final_for-analysis.txt", quote = F, sep = "\t", row.names = T)

#------------------------------------------------#
# subset by mile (i.e. elevation band/community) #
#------------------------------------------------#

final <- read.table("final_for-analysis.txt", sep = "\t", header = T)

sub_m1 <- subset(final, Mile == 1);
sub_m2 <- subset(final, Mile == 2);
sub_m3 <- subset(final, Mile == 3);
sub_m4 <- subset(final, Mile == 4);
sub_m5 <- subset(final, Mile == 5);


write.table(sub_m1, file = "sub_m1.txt", quote = F, sep = "\t", row.names = T)
write.table(sub_m2, file = "sub_m2.txt", quote = F, sep = "\t", row.names = T)
write.table(sub_m3, file = "sub_m3.txt", quote = F, sep = "\t", row.names = T)
write.table(sub_m4, file = "sub_m4.txt", quote = F, sep = "\t", row.names = T)
write.table(sub_m5, file = "sub_m5.txt", quote = F, sep = "\t", row.names = T)

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


## Sources overlap functions ##
source(".../overlap.R")


#------------------------------------------------------------------------------#
# apply overlap function to dataset subset per elevation band (i.e. community) #
#------------------------------------------------------------------------------#

## mile 1 ##
pair_overlap_m1 <- as.data.frame(t(overlap2(mydata = "sub_m1.txt"))) #data.frame with pairs of species in line and overlap coefficients in column

pair_overlap_m1 <- rownames_to_column(pair_overlap_m1, var = "sp_year") #add name to row.names column

overlap_m1 <- pair_overlap_m1 %>%
              separate(sp_year, c("pair", "year"), "-") #split columns with pair of species and year

overlap_m1['mile']='1' #add column "mile"

write.table(overlap_m1, file = "overlap_tot_m1.txt", quote = F, sep = "\t", row.names = T)

overlap_filt1 <- overlap_m1 %>% group_by(pair) %>% filter(n() >= 10 & length(year) >= 20) #filter pairs of species that occur at least 10 years over 20 a year period

overlap_filt1$year <- as.numeric(overlap_filt1$year)

overlap_filt1$centered_year <- scale(overlap_filt1$year, scale = FALSE) #add column with centered years
overlap_filt1$centered_year_sq <-as.numeric((overlap_filt1$centered_year)^2) #add column with squared centered years

overlap_filt1 <- overlap_filt1[c("pair", "mile", "year", "centered_year", "centered_year_sq", "overlap_Est")] #reorder columns

write.table(overlap_filt1, file = "overlap_filt_m1.txt", quote = F, sep = "\t", row.names = T)


## mile 2 ##
pair_overlap_m2 <- as.data.frame(t(overlap2(mydata = "sub_m2.txt"))) #data.frame with pairs of species in line and overlap coefficients in column

pair_overlap_m2 <- rownames_to_column(pair_overlap_m2, var = "sp_year") #add name to row.names column

overlap_m2 <- pair_overlap_m2 %>%
  separate(sp_year, c("pair", "year"), "-") #split columns with pair of species and year

overlap_m2['mile']='2' #add column "mile"

write.table(overlap_m2, file = "overlap_tot_m2.txt", quote = F, sep = "\t", row.names = T)

overlap_filt2 <- overlap_m2 %>% group_by(pair) %>% filter(n() >= 10 & length(year) >= 20) #filter pairs of species that occur at least 10 years over 20 a year period

overlap_filt2$year <- as.numeric(overlap_filt2$year)

overlap_filt2$centered_year <- scale(overlap_filt2$year, scale = FALSE) #add column with centered years
overlap_filt2$centered_year_sq <-as.numeric((overlap_filt2$centered_year)^2) #add column with squared centered years

overlap_filt2 <- overlap_filt2[c("pair", "mile", "year", "centered_year", "centered_year_sq", "overlap_Est")] #reorder columns

write.table(overlap_filt2, file = "overlap_filt_m2.txt", quote = F, sep = "\t", row.names = T)


## mile 3 ##
pair_overlap_m3 <- as.data.frame(t(overlap2(mydata = "sub_m3.txt"))) #data.frame with pairs of species in line and overlap coefficients in column

pair_overlap_m3 <- rownames_to_column(pair_overlap_m3, var = "sp_year") #add name to row.names column

overlap_m3 <- pair_overlap_m3 %>%
  separate(sp_year, c("pair", "year"), "-") #split columns with pair of species and year

overlap_m3['mile']='3' #add column "mile"

write.table(overlap_m3, file = "overlap_tot_m3.txt", quote = F, sep = "\t", row.names = T)

overlap_filt3 <- overlap_m3 %>% group_by(pair) %>% filter(n() >= 10 & length(year) >= 20) #filter pairs of species that occur at least 10 years over 20 a year period

overlap_filt3$year <- as.numeric(overlap_filt3$year)

overlap_filt3$centered_year <- scale(overlap_filt3$year, scale = FALSE) #add column with centered years
overlap_filt3$centered_year_sq <-as.numeric((overlap_filt3$centered_year)^2) #add column with squared centered years

overlap_filt3 <- overlap_filt3[c("pair", "mile", "year", "centered_year", "centered_year_sq", "overlap_Est")] #reorder columns

write.table(overlap_filt3, file = "overlap_filt_m3.txt", quote = F, sep = "\t", row.names = T)


## mile 4 ##
pair_overlap_m4 <- as.data.frame(t(overlap2(mydata = "sub_m4.txt"))) #data.frame with pairs of species in line and overlap coefficients in column

pair_overlap_m4 <- rownames_to_column(pair_overlap_m4, var = "sp_year") #add name to row.names column

overlap_m4 <- pair_overlap_m4 %>%
  separate(sp_year, c("pair", "year"), "-") #split columns with pair of species and year

overlap_m4['mile']='4' #add column "mile"

write.table(overlap_m4, file = "overlap_tot_m4.txt", quote = F, sep = "\t", row.names = T)

overlap_filt4 <- overlap_m4 %>% group_by(pair) %>% filter(n() >= 10 & length(year) >= 20) #filter pairs of species that occur at least 10 years over 20 a year period

overlap_filt4$year <- as.numeric(overlap_filt4$year)

overlap_filt4$centered_year <- scale(overlap_filt4$year, scale = FALSE) #add column with centered years
overlap_filt4$centered_year_sq <-as.numeric((overlap_filt4$centered_year)^2) #add column with squared centered years

overlap_filt4 <- overlap_filt4[c("pair", "mile", "year", "centered_year", "centered_year_sq", "overlap_Est")] #reorder columns

write.table(overlap_filt4, file = "overlap_filt_m4.txt", quote = F, sep = "\t", row.names = T)


## mile 5 ##
pair_overlap_m5 <- as.data.frame(t(overlap2(mydata = "sub_m5.txt"))) #data.frame with pairs of species in line and overlap coefficients in column

pair_overlap_m5 <- rownames_to_column(pair_overlap_m5, var = "sp_year") #add name to row.names column

overlap_m5 <- pair_overlap_m5 %>%
  separate(sp_year, c("pair", "year"), "-") #split columns with pair of species and year

overlap_m5['mile']='5' #add column "mile"

write.table(overlap_m5, file = "overlap_tot_m5.txt", quote = F, sep = "\t", row.names = T)

overlap_filt5 <- overlap_m5 %>% group_by(pair) %>% filter(n() >= 10 & length(year) >= 20) #filter pairs of species that occur at least 10 years over 20 a year period

overlap_filt5$year <- as.numeric(overlap_filt5$year)

overlap_filt5$centered_year <- scale(overlap_filt5$year, scale = FALSE) #add column with centered years
overlap_filt5$centered_year_sq <-as.numeric((overlap_filt5$centered_year)^2) #add column with squared centered years

overlap_filt5 <- overlap_filt5[c("pair", "mile", "year", "centered_year", "centered_year_sq", "overlap_Est")] #reorder columns

write.table(overlap_filt5, file = "overlap_filt_m5.txt", quote = F, sep = "\t", row.names = T)



##bind the 5 miles together
overlap_filt_tot <- rbind(overlap_filt1, overlap_filt2, overlap_filt3, overlap_filt4, overlap_filt5)

write.table(overlap_filt_tot, file = "overlap_all_miles.txt", quote = F, sep = "\t", row.names = T)



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

# ----------------------------------------------------------------------------------#
# Function to estimate flowering overlap between elevation bands (i.e. communities) #
# ----------------------------------------------------------------------------------#

library(dummy)
library(igraph)
library(tidyr)
library(dplyr)
library(overlap)

## function to obtain a graph + matrix
{
my.graph.fun<-function(el){
  mat<-as.matrix(el,dim(el)[1],2)
  mat.mod<-as.matrix(cbind(sapply(mat[,1],FUN = function(x) paste0("M_",x)),sapply(mat[,2],FUN = function(x) paste0("P_",x))))
  graph_from_edgelist(mat.mod,directed=F)
}

## function to obtain an incidence matrix
my.incidence<-function(gr){
  adj<-as.data.frame(as.matrix(as_adjacency_matrix(gr)))
  adj<-adj[order(names(adj)),order(names(adj))]
  M.index<-max(grep("M_",names(adj)))
  adj[1:M.index,(M.index+1):length(names(adj))]
}
}

#function to subset each species present in both miles in each year
#mydata<-read.table("m12.txt",sep="\t",header=T)
#mydata<-read.table("m13.txt",sep="\t",header=T)
#mydata<-read.table("m14.txt",sep="\t",header=T)
#mydata<-read.table("m15.txt",sep="\t",header=T)
#mydata<-read.table("m23.txt",sep="\t",header=T)
#mydata<-read.table("m24.txt",sep="\t",header=T)
#mydata<-read.table("m25.txt",sep="\t",header=T)
#mydata<-read.table("m34.txt",sep="\t",header=T)
#mydata<-read.table("m35.txt",sep="\t",header=T)
mydata<-read.table("m45.txt",sep="\t",header=T)


filtersp<-function(mydata)
{
  mydata<-mydata
  #select unique years
  y<-unique(mydata$Year)
  #for each year
  for(annee in y)
  {
    #print("annee")
    #print(annee)
    #select all lines in the table 'mydata' that correspond to the year annee
    subset_year<-mydata[which(mydata$Year == annee),]
    #create the filename
    filename<-subset_year$Year[1]
    filenameTXT<-paste(filename, "_Y.txt", sep="")
#mydata<-subset(mydata, Year==1984)
    subsp<-subset_year[,1:2]
    subsp<-my.incidence(my.graph.fun(subsp))

#binary all species
    subsp<-t(subsp)
    subsp <- subsp
    subsp[subsp>=1]<-1

    subsp <- cbind(subsp, Total = rowSums(subsp[,1:2])) #add sum column
    sp.names<- rownames(subsp) #save names as character
    rownames(subsp) <- NULL #remove rownames 
    subsp <-subsp[,c(3,1,2)] #1st column Tot
    data0 <- as.data.frame(cbind(subsp,sp.names)) #use rownames as column
    data0<-data0[do.call(order, c(data0, list(decreasing=TRUE))),] #decrease order, all the sum = 0 are at the end of the table

    subsp<-data0[which(data0$Total == 2),]

#mydata.sp<-data0[1:203,] #keep just the lines with tot =2
    subsp<-subsp[,c(4,1:3)]
    subsp<-subsp[do.call(order, c(subsp, list(decreasing=FALSE))),] #decrease order, all the sum = 0 are at the end of the table

goodsp<-subsp %>% separate(sp.names,c("type","species"),"_")
goodsp<-goodsp$species

myfilter<-subset_year %>% filter(Species %in% goodsp)
write.table(myfilter,file = filenameTXT,quote = F,sep = "\t",row.names = T)

  }
}
df <- filtersp(mydata = mydata)

#Load all the txt per Y and create one dataset (species repeated x33 years)
{
  filelist = list.files(pattern = "Y.txt")
  datalist = lapply(filelist, function(x)read.table(x, header=T,sep="\t"))
  datafr = do.call("rbind", datalist) 
  write.table(datafr,file = "datafr_subsp.txt",quote = F,sep = "\t")
}

#Function overlap subset year + subset same species in different miles
overlap3<-function(mydata)
{
  results<-NULL
  mydata<-mydata
  #select unique years
  y<-unique(mydata$Year)
  #for each year
  for(annee in y)
  {
    #print("annee")
    #print(annee)
    #select all lines in the table 'mydata' that correspond to the year annee
    subset_year<-mydata[which(mydata$Year == annee),]
    sp<-unique(subset_year$Species)#select species from the year "annee" (of the loop)
    espece_done<-NULL #vector that contains species already compared (to avoid replicates). It is empty, it will fill as comparisons occur
    m<-unique(subset_year$Mile) #select miles from the year "annee" (of the loop)
    #for every espece1
    for(espece in sp)
    {
      espece_done<-cbind(espece_done, espece)
      #select all lines of the sub-table "annee" that correspond to the species specie
      #select the species to be compared with
      subset_specie<- subset_year[which(subset_year$Species == espece),] #subset of espece1 within a year
      #____#
      # se subset_specie$mile ha sia m1 che m2 allora continua con la funzione di seguito
      #altrimenti passa alla specie successiva
      #____#
      for(miglia in m)
      {
        m1<- subset_specie[which(subset_specie$Mile == miglia),]
        m2 <- subset_specie[- which(subset_specie$Mile == miglia),]
        if((m1$Species) == (m2$Species))
        {
          if(length(m1$Species)[1] > 50  & length(m2$Species)[1] > 50)
          {
            #print("risotto")
            res1<-as.matrix(overlapEst( m1$doyRad,  m2$doyRad, type = "Dhat4")) #result of the function overlapEst
            colnames(res1)<-paste(paste(espece, espece, sep="_and_"), annee, sep="_*_") #add column name with different separators
            results<-cbind(results, res1) #add this result to the final table
          }
          else
          {
            #print("pizza")
            res2<-as.matrix(overlapEst( m1$doyRad,  m2$doyRad, type = "Dhat1")) #result of the function overlapEst
            colnames(res2)<-paste(paste(espece, espece, sep="_and_"), annee, sep="_*_") #add column name with different separators
            results<-cbind(results, res2) #add this result to the final tableresults<-cbind(results, res) #add this result to the final table
          }
        }
      }
    }
  }
  return(results)
}

mydata<-read.table("datafr_subsp.txt",sep="\t",header=T)
df<-overlap3(mydata = mydata)
df<-as.data.frame(t(df))
df$double<-rep(1:2)
df1<-df[-which(df$double == 2),]
species<- rownames(df1) #save names as character
rownames(df1) <- NULL
df1$species<-species
df1<-df1[,c(3,1)]


#add columns with centered year and centered year squared
df1 <- df1 %>%
  separate(species, c("pair", "year"), "_._") #split columns with pair of species and year
df1$year <- as.numeric(df1$year)
df1$center_scale_year <- scale(df1$year, scale = TRUE) #add column with centered years
df1$center_scale_year_sq <-as.numeric((df1$center_scale_year)^2) #add column with squared centered years
df1['mile_pair']='45' #add column with "mile_pair"
names(df1)[names(df1)=="Dhat1"] <- "overlap_Est" #rename overlap column
df1 <- df1[c("pair", "mile_pair", "year", "center_scale_year", "center_scale_year_sq", "overlap_Est")] #reorder columns



#write.table(df1,file = "m12_overlap.txt",quote = F,sep = "\t",row.names = T)
#write.table(df1,file = "m13_overlap.txt",quote = F,sep = "\t",row.names = T)
#write.table(df1,file = "m14_overlap.txt",quote = F,sep = "\t",row.names = T)
#write.table(df1,file = "m15_overlap.txt",quote = F,sep = "\t",row.names = T)
#write.table(df1,file = "m23_overlap.txt",quote = F,sep = "\t",row.names = T)
#write.table(df1,file = "m24_overlap.txt",quote = F,sep = "\t",row.names = T)
#write.table(df1,file = "m25_overlap.txt",quote = F,sep = "\t",row.names = T)
#write.table(df1,file = "m34_overlap.txt",quote = F,sep = "\t",row.names = T)
#write.table(df1,file = "m35_overlap.txt",quote = F,sep = "\t",row.names = T)
write.table(df1,file = "m45_overlap.txt",quote = F,sep = "\t",row.names = T)

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

#------------------------------------------------------#
# Plot flowering phenology overlap between two species #
#------------------------------------------------------#

library(overlap)

m1<-read.table("sub_m1.txt",header=T,sep="\t")

m1_sub <- subset(m1, Year == 1989, select = c(Species, doyRad)) #Subset year 1989
table(m1_sub$Species)

# Example of high flowering overlap
heri <- subset(m1_sub, Species =='Herissantia crispa')
plot1 <- densityPlot(heri$doyRad, xscale = 365, rug=TRUE, adjust =1)

just <- subset(m1_sub, Species =='Justicia longii')
plot2 <- densityPlot(just$doyRad, xscale = 365, rug=TRUE, adjust =1)

min(length(heri), length(just))

heri_just <- overlapEst(heri$doyRad, just$doyRad, type="Dhat1")
heri_just

plot_overlap <- overlapPlot(heri$doyRad, just$doyRad, xscale = 365, linewidth = c(2, 2), xaxt = "n", yaxt = "n", rug = TRUE, axes = F, cex.lab = 1.5, main="", ylab = "Density of flowering", xlab="Day of year")
box(bty="l")
axis(1,cex.axis = 1.5)
axis(2,cex.axis = 1.5)
legend(-30, 0.01, c("Herissantia crispa", "Justicia longii"), lty = c(1,2), col = c("black", "blue"), lwd = 2, pt.cex = 1, cex = 2, bty='n', y.intersp = 0.1)



# Example of low flowering overlap
carn <- subset(m1_sub, Species =='Carnegiea gigantea')
plot3 <- densityPlot(carn$doyRad, xscale = 365, rug=TRUE, adjust =1)

medi <- subset(m1_sub, Species =='Medicago polymorpha')
plot4 <- densityPlot(medi$doyRad, xscale = 365, rug=TRUE, adjust =1)


carn_medi <- overlapEst(carn$doyRad, medi$doyRad, type="Dhat1")
carn_medi

plot_overlap <- overlapPlot(carn$doyRad, medi$doyRad, xscale = 365, linewidth = c(2, 2), xaxt = "n", yaxt = "n", rug = TRUE, axes = F, cex.lab = 1.5, main="", ylab = "Density of flowering", xlab="Day of year")
box(bty="l")
axis(1,cex.axis = 1.5)
axis(2,cex.axis = 1.5)
legend(120, 0.018, c("Carnegiea gigantea", "Medicago polymorpha"), lty = c(1,2), col = c("black", "blue"), lwd = 2, pt.cex = 1, cex = 2, bty='n', y.intersp = 0.1)


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

#---------------------------------------------------------------------------------------#
# Zero-and-one inflated beta regressions using Bayesian generalized linear mixed models # 
#---------------------------------------------------------------------------------------#

library('brms') #load package


## Within elevation band (i.e. community)


model_int <- function()  # global model to test for interactions between elevation bands (i.e. community) and year
{
  inter.zoib <- brm(brmsformula(overlap_Est ~ center_scale_year + mile + center_scale_year:mile + (1|pair), phi ~ mile, zoi ~ mile, coi ~ mile, family=zero_one_inflated_beta()), data=overlap_all_miles, iter=10000, cores = 4)
  
  saveRDS(inter.zoib, "inter.zoib.RDS")
  
}
model_int()



models_mile1 <- function()  # model to test for effect of year and year squared within elevation band (i.e. community). Use same model for each elevation band
{
  m1y.zoib <- brm(brmsformula(overlap_Est ~ center_scale_year + (1|pair), phi ~ center_scale_year + (1|pair), zoi ~ center_scale_year + (1|pair), coi ~center_scale_year + (1|pair), family=zero_one_inflated_beta()), data=overlap_m1, iter=10000, cores = 4)
  #m1y.zoib_sq <- brm(brmsformula(overlap_Est ~ center_scale_year_sq + (1|pair), phi ~ center_scale_year_sq + (1|pair), zoi ~ center_scale_year_sq + (1|pair), coi ~ center_scale_year_sq + (1|pair), family=zero_one_inflated_beta()), data=overlap_m1, iter=10000, cores = 4)
    
  saveRDS(m1y.zoib, "m1y.zoib.RDS")
  #saveRDS(m1y.zoib_sq, "m1y.zoib_sq.RDS")
    
}
models_mile1()


# Model comparison
loo_1 <- loo(m1y.zoib_sq, m1y.zoib, nsamples = 2000) #set nsamples lower if too computationally demanding
loo_2 <- loo(m2y.zoib_sq, m2y.zoib, nsamples = 2000)
loo_3 <- loo(m3y.zoib_sq, m3y.zoib, nsamples = 2000)
loo_4 <- loo(m4y.zoib_sq, m4y.zoib, nsamples = 2000)
loo_5 <- loo(m5y.zoib_sq, m5y.zoib, nsamples = 2000)


# Model summary and checks. Use same model for each elevation band
m1<-readRDS("m1y.zoib.RDS")
summary(m1)
print(summary(m1), digits = 4)

cef1 <- conditional_effects(m1, dpar = "mu")
cef1_year <- cef1$center_scale_year[,5] #retrieve effect (year)
cef1_pred <- cef1$center_scale_year[,6] #retrieve predicted values of the response
cef1_lowCI <- cef1$center_scale_year[,8] #retrieve lower CI
cef1_highCI <- cef1$center_scale_year[,9] #retrieve higher CI
m1_pred<-data.frame(cef1_year, cef1_pred, cef1_lowCI, cef1_highCI)
write.table(m1_pred, file = "m1_pred.txt", quote = F, sep = "\t")

#plots
trace.dens1 <- plot(m1) #Trace and Density Plots for MCMC Samples
ppcheck1 <- pp_check(m1) #Posterior predictive checks
qq.plot1 <- pp_check(m1, type = "loo_pit") #Q-Q plot
pred_m1 <- plot(conditional_effects(m1, dpar = "mu"), #Predicted values
                points = F, 
                point_args = list(width = .05, shape = 1))




## Between elevation bands (i.e. communities)

model_int <- function() # global model to test for interactions between elevation band (i.e. community) pairs and year
{
  inter.btwn.zoib <- brm(brmsformula(overlap_Est ~ center_scale_year + mile_pair + center_scale_year:mile_pair + (1|pair), phi ~ mile_pair, zoi ~ mile_pair, coi ~ mile_pair, family=zero_one_inflated_beta()), data=btwn_miles, iter=10000, cores = 4)
  
  saveRDS(inter.btwn.zoib, "inter.btwn.zoib.RDS")
  
}
model_int()



models_miles15 <- function() # models to test for effect of year and year squared  between elevation bands (i.e. communities) 1 through 5. Use same model for each elevation band pair.
{
  m12y.zoib <- brm(brmsformula(overlap_Est ~ center_scale_year + (1|pair), phi ~ center_scale_year, zoi ~ center_scale_year, coi ~center_scale_year, family=zero_one_inflated_beta()), data=m12, iter=10000, cores = 4)
  m12y.zoib_sq <- brm(brmsformula(overlap_Est ~ center_scale_year_sq + (1|pair), phi ~ center_scale_year_sq, zoi ~ center_scale_year_sq, coi ~ center_scale_year_sq, family=zero_one_inflated_beta()), data=m12, iter=10000, cores = 4)
  
  m13y.zoib <- brm(brmsformula(overlap_Est ~ center_scale_year + (1|pair), phi ~ center_scale_year, zoi ~ center_scale_year, coi ~center_scale_year, family=zero_one_inflated_beta()), data=m13, iter=10000, cores = 4)
  m13y.zoib_sq <- brm(brmsformula(overlap_Est ~ center_scale_year_sq + (1|pair), phi ~ center_scale_year_sq, zoi ~ center_scale_year_sq, coi ~ center_scale_year_sq, family=zero_one_inflated_beta()), data=m13, iter=10000, cores = 4)
  
  m14y.zoib <- brm(brmsformula(overlap_Est ~ center_scale_year + (1|pair), phi ~ center_scale_year, zoi ~ center_scale_year, coi ~center_scale_year, family=zero_one_inflated_beta()), data=m14, iter=10000, cores = 4)
  m14y.zoib_sq <- brm(brmsformula(overlap_Est ~ center_scale_year_sq + (1|pair), phi ~ center_scale_year_sq, zoi ~ center_scale_year_sq, coi ~ center_scale_year_sq, family=zero_one_inflated_beta()), data=m14, iter=10000, cores = 4)
  
  m15y.zoib <- brm(brmsformula(overlap_Est ~ center_scale_year + (1|pair), phi ~ center_scale_year, zoi ~ center_scale_year, coi ~center_scale_year, family=zero_one_inflated_beta()), data=m15, iter=10000, cores = 4)
  m15y.zoib_sq <- brm(brmsformula(overlap_Est ~ center_scale_year_sq + (1|pair), phi ~ center_scale_year_sq, zoi ~ center_scale_year_sq, coi ~ center_scale_year_sq, family=zero_one_inflated_beta()), data=m15, iter=10000, cores = 4)
  
  saveRDS(m12y.zoib, "m12y.zoib.RDS")
  saveRDS(m12y.zoib_sq, "m12y.zoib_sq.RDS")
  
  saveRDS(m13y.zoib, "m13y.zoib.RDS")
  saveRDS(m13y.zoib_sq, "m13y.zoib_sq.RDS")
  
  saveRDS(m14y.zoib, "m14y.zoib.RDS")
  saveRDS(m14y.zoib_sq, "m14y.zoib_sq.RDS")
  
  saveRDS(m15y.zoib, "m15y.zoib.RDS")
  saveRDS(m15y.zoib_sq, "m15y.zoib_sq.RDS")
}
models_miles15()


# Model comparison
loo_12 <- loo(m12y.zoib_sq, m12y.zoib, nsamples = 2000) #set nsamples if too computationally demanding
loo_13 <- loo(m13y.zoib_sq, m13y.zoib, nsamples = 2000)
loo_14 <- loo(m14y.zoib_sq, m14y.zoib, nsamples = 2000)
loo_15 <- loo(m15y.zoib_sq, m15y.zoib, nsamples = 2000)
loo_23 <- loo(m23y.zoib_sq, m23y.zoib, nsamples = 2000)
loo_24 <- loo(m24y.zoib_sq, m24y.zoib, nsamples = 2000)
loo_25 <- loo(m25y.zoib_sq, m25y.zoib, nsamples = 2000)
loo_34 <- loo(m34y.zoib_sq, m34y.zoib, nsamples = 2000)
loo_35 <- loo(m35y.zoib_sq, m35y.zoib, nsamples = 2000)
loo_45 <- loo(m45y.zoib_sq, m45y.zoib, nsamples = 2000)

# Model summary and checks. Use same model for each elevation band pair
m12<-readRDS("m12y.zoib_sq.RDS")
summary(m12)
print(summary(m12), digits = 4)
R2_12 <- loo_R2(m12) #R-squared

cef12 <- conditional_effects(m12, dpar = "mu")
cef12_year <- cef12$center_scale_year[,5] #retrieve effect (year)
cef12_pred <- cef12$center_scale_year[,6] #retrieve predicted values of the response
cef12_lowCI <- cef12$center_scale_year[,8] #retrieve lower CI
cef12_highCI <- cef12$center_scale_year[,9] #retrieve higher CI
m12_pred<-data.frame(cef12_year, cef12_pred, cef12_lowCI, cef12_highCI)
write.table(m12_pred, file = "m12_pred.txt", quote = F, sep = "\t")

#plots
trace.dens12 <- plot(m12) #Trace and Density Plots for MCMC Samples
ppcheck12 <- pp_check(m12) #posterior predictive checks
qq.plot12 <- pp_check(m12, type = "loo_pit", nsamples = 2000) #Q-Q plot
pred_m12 <- plot(conditional_effects(m12, dpar = "mu"), #predicted values
                points = F, 
                point_args = list(width = .05, shape = 1))



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



#---------------------------------------------------------------------------------#
# Models to estimate effect of climate on flowering overlap within elevation band #
#---------------------------------------------------------------------------------#


model_int <- function() #Test for effect of temperature and precipitation on overlap. Use for each elevation band
{
  overlap_climate_m1 <- brm(brmsformula(overlap_Est ~ center_scale_temp + center_scale_prec + center_scale_temp:center_scale_prec + (1|pair), phi ~ center_scale_temp*center_scale_prec + (1|pair), zoi ~ center_scale_temp*center_scale_prec + (1|pair), coi ~ center_scale_temp*center_scale_prec + (1|pair), family=zero_one_inflated_beta()), data = subset(overlap_climate, mile == 1), iter = 10000, cores = 4)
  
  saveRDS(overlap_climate_m1, "overlap.climate.m1b.zoib.RDS")
  
}
model_int()





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

#---------------------------------------------------------------------------------------------------------------------------#
# Models to test for temporal changes of climatic variables (temperature and precipitation) across the elevational gradient #
#---------------------------------------------------------------------------------------------------------------------------#

library(dplyr)
library(ggplot2)

rain <- read.table("prism_precipitation.txt", sep = "\t", header = T)
temp <- read.table("prism_temperature.txt", sep = "\t", header = T)


## PRECIPITATION ##

##Calculate summary values
sem <- function(x) sd(x)/sqrt(length(x)) #function to calculate standard error

rain_avg <- group_by(rain, year, elevation) %>%
                 summarise(
                  count = n(),
                  sum = sum(ppt_mm_gauge, na.rm = TRUE),
                  mean = mean(ppt_mm_gauge, na.rm = TRUE),
                  sd = sd(ppt_mm_gauge, na.rm = TRUE),
                  sem = sem(ppt_mm_gauge),
                  median = median(ppt_mm_gauge, na.rm = TRUE)
              )

write.table(rain_avg, "precip_summary.txt", sep = "\t")

##Plot precipitation by elevation and year
line_rain <- ggplot(rain_avg, aes(x = year, y = sum, color = factor(elevation))) +
  #geom_line(aes(color = factor(elevation)), size = 1.2) +
  geom_point(size = 2) +
  geom_smooth(method='lm', alpha = .15, size = 2, aes(fill = -elevation)) +
  scale_colour_manual(values = c("skyblue1", "deepskyblue", "deepskyblue4")) +
  scale_x_continuous(breaks=seq(1984,2019,5)) +
  labs(x = "\nYear", y = expression(paste("Precipitation", " ", "(mm year" ^-{1}, ")"))) +
  theme(panel.background = element_blank(), axis.line.x = element_line(color="black", size = 0.5),axis.line.y = element_line(color="black", size = 0.5)) +
  theme(axis.text=element_text(size=30), axis.title = element_text(size=30)) +
  theme(legend.position = "none")

line_rain

##Linear models
library(rcompanion)

rain_sum <- read.table("precip_summary.txt", sep = "\t", header = T)
rain_sum$elevation <- as.factor(rain_sum$elevation)
rain_sum$log_sum <- log(rain_sum$sum)

plotNormalHistogram(rain_sum$sum)
plotNormalHistogram(rain_sum$log_sum)

lm1 <- lm(sum ~ year * elevation, data = rain_sum)
anova(lm1)
summary(lm1)
plotNormalHistogram(resid(lm1))
qqnorm(resid(lm1))
qqline(resid(lm1))

lm2 <- lm(log_sum ~ year * elevation, data = rain_sum) #best fit
anova(lm2)
summary(lm2)
plotNormalHistogram(resid(lm2))
qqnorm(resid(lm2))
qqline(resid(lm2))

#because interaction not significant in the previous lm, lmm with elevation as random effect to allow intercepts to vary
library(lme4)
library(lmerTest)
library(MuMIn)
lme1 <- lmer(log_sum ~ year * elevation + (1 | elevation), data = rain_sum) 
anova(lme1)
summary(lme1)
r.squaredGLMM(lme1)


## TEMPERATURE ##

##Calculate summary values
sem <- function(x) sd(x)/sqrt(length(x)) #function to calculate standard error

temp_avg <- group_by(temp, year, elevation) %>%
  summarise(
    count = n(),
    mean = mean(tmp, na.rm = TRUE),
    sd = sd(tmp, na.rm = TRUE),
    sem = sem(tmp),
    median = median(tmp, na.rm = TRUE)
  )

write.table(temp_avg, "temp_summary.txt", sep = "\t")

##Plot temperature by elevation and year
temp_avg <- read.table("temp_summary.txt", sep = "\t", header = T)
temp_avg$elevation <- as.factor(temp_avg$elevation)

line_temp <- ggplot(temp_avg, aes(x = year, y = mean, color = elevation)) +
  geom_smooth(method='lm', alpha = .15, size = 2, aes(fill = elevation)) +
  geom_point(size = 2) +
  scale_colour_manual(values = c("red", "orange", "gold")) +
  scale_x_continuous(breaks=seq(1984,2019,5)) +
  labs(x = "\nYear", y = "Temperature °C\n") +
  theme(panel.background = element_blank(), axis.line.x = element_line(color="black", size = 0.5),axis.line.y = element_line(color="black", size = 0.5)) +
  theme(axis.text=element_text(size=30), axis.title = element_text(size=30)) +
  theme(legend.position = "none")

line_temp

##Linear models
library(rcompanion)

temp_avg <- read.table("temp_summary.txt", sep = "\t", header = T)
temp_avg$elevation <- as.factor(temp_avg$elevation)
temp_avg$log_mean <- log(temp_avg$mean)

plotNormalHistogram(temp_avg$mean)
plotNormalHistogram(temp_avg$log_mean)

lm3 <- lm(mean ~ year * elevation, data = temp_avg)
anova(lm3)
summary(lm3)
plotNormalHistogram(resid(lm3))
qqnorm(resid(lm3))
qqline(resid(lm3))

lm4 <- lm(log_mean ~ year * elevation, data = temp_avg) #best fit
anova(lm4)
summary(lm4)
plotNormalHistogram(resid(lm4))
qqnorm(resid(lm4))
qqline(resid(lm4))

lm5 <- lm(log_mean ~ year, data = subset(temp_avg, elevation == 1)) #elevation 1
anova(lm5)
summary(lm5)
plotNormalHistogram(resid(lm5))

lm6 <- lm(log_mean ~ year, data = subset(temp_avg, elevation == 2)) #elevation 2
anova(lm6)
summary(lm6)
plotNormalHistogram(resid(lm6))

lm7 <- lm(log_mean ~ year, data = subset(temp_avg, elevation == 3)) #elevation 3
anova(lm7)
summary(lm7)
plotNormalHistogram(resid(lm7))