Supplementary file Code S3 - Multiple spatial and temporal replicates 

# Community dynamics analysis of the BCI dataset
# Data for the first eight censuses were downloaded. Only alive, main stems were downloaded. 
# Users must request the data from: http://ctfs.si.edu/webatlas/datasets/bci/
# Further description of this worked case study example is provided in Supplemental Article S2.

# Note that the code for the randomisations to test for significance takes quite a long time.

# Install required packages
library(tidyverse) # for data handling and manipulation
library(vegan) # for a range of community analyses
library(adespatial)
library(spatstat)
library(codyn) # to calculate species turnover
library(zetadiv) # to calculate zeta diversity

# Read in the datasets and add a column for most records were taken
raw_census1 <- read.table("PlotDataReport04-22-2018_census-1.txt", header = T, sep = "\t")
  raw_census1$Year <- 1982
  raw_census1 <- raw_census1 %>% filter(!is.na(PX)) %>% filter(Latin != "Unidentified species")
raw_census2 <- read.table("PlotDataReport04-22-2018_census-2.txt", header = T, sep = "\t")
  raw_census2$Year <- 1985
  raw_census2 <- raw_census2 %>% filter(!is.na(PX)) %>% filter(Latin != "Unidentified species")
raw_census3 <- read.table("PlotDataReport04-22-2018_census-3.txt", header = T, sep = "\t")
  raw_census3$Year <- 1990
  raw_census3 <- raw_census3 %>% filter(!is.na(PX)) %>% filter(Latin != "Unidentified species")
raw_census4 <- read.table("PlotDataReport04-22-2018_census-4.txt", header = T, sep = "\t")
  raw_census4$Year <- 1995
  raw_census4 <- raw_census4 %>% filter(!is.na(PX)) %>% filter(Latin != "Unidentified species")
raw_census5 <- read.table("PlotDataReport04-22-2018_census-5.txt", header = T, sep = "\t")
  raw_census5$Year <- 2000
  raw_census5 <- raw_census5 %>% filter(!is.na(PX)) %>% filter(Latin != "Unidentified species")
raw_census6 <- read.table("PlotDataReport04-22-2018_census-6.txt", header = T, sep = "\t")
  raw_census6$Year <- 2005
  raw_census6 <- raw_census6 %>% filter(!is.na(PX)) %>% filter(Latin != "Unidentified species")
raw_census7 <- read.table("PlotDataReport04-22-2018_census-7.txt", header = T, sep = "\t")
  raw_census7$Year <- 2010
  raw_census7 <- raw_census7 %>% filter(!is.na(PX)) %>% filter(Latin != "Unidentified species")
raw_census8 <- read.table("PlotDataReport04-22-2018_census-8.txt", header = T, sep = "\t")
  raw_census8$Year <- 2015
  raw_census8 <- raw_census8 %>% filter(!is.na(PX)) %>% filter(Latin != "Unidentified species")

# Merge datasets and sort dataframe
bci <- rbind(raw_census1, raw_census2, raw_census3, raw_census4, raw_census5, raw_census6, raw_census7, raw_census8)
bci$presence <- 1

# Remove raw data
#rm(raw_census1, raw_census2, raw_census3, raw_census4, raw_census5, raw_census6, raw_census7, raw_census8)

##########################################################
# Exploratory Data Analysis (EDA)
##########################################################

# Count the number of species and samples
length(unique(bci$Latin)) # eight censuses and 324 species
unique(bci$Latin) # list the 324 species

# Calculate the total abundance of each species across all censuses  
bci_census_summary <- bci %>% group_by(Year, Latin) %>% summarise(Abundance = sum(presence)) 
bci_census_summary %>% ggplot(aes(x = Year, y = Abundance, col = Latin)) + geom_line(show.legend = F) + theme_bw()
bci_census_summary %>% ggplot(aes(x = Abundance)) + geom_histogram(bins = 25) + theme_bw()

##########################################################
# Prepare datasets at different spatial scales
##########################################################

# census years
census_years <- c(1982, 1985, 1990, 1995, 2000, 2005, 2010, 2015)

# set plot coordinates
xmin=0; xmax=1000; ymin=0; ymax=500  # plot dimensions

# quadrat sizes to generate datasets for
quadrat_sizes <- c(20, 50, 125, 250, 500)

# ppp to quadrat-level community datasets
for(k in 1:8){ # census loop
  if(k == 1) { ppp.dat <- ppp(x = raw_census1$PX, y = raw_census1$PY, marks = raw_census1$Latin, window = owin(c(0, 1000), c(0, 500))) }
  if(k == 2) { ppp.dat <- ppp(x = raw_census2$PX, y = raw_census2$PY, marks = raw_census2$Latin, window = owin(c(0, 1000), c(0, 500))) }
  if(k == 3) { ppp.dat <- ppp(x = raw_census3$PX, y = raw_census3$PY, marks = raw_census3$Latin, window = owin(c(0, 1000), c(0, 500))) }
  if(k == 4) { ppp.dat <- ppp(x = raw_census4$PX, y = raw_census4$PY, marks = raw_census4$Latin, window = owin(c(0, 1000), c(0, 500))) }
  if(k == 5) { ppp.dat <- ppp(x = raw_census5$PX, y = raw_census5$PY, marks = raw_census5$Latin, window = owin(c(0, 1000), c(0, 500))) }
  if(k == 6) { ppp.dat <- ppp(x = raw_census6$PX, y = raw_census6$PY, marks = raw_census6$Latin, window = owin(c(0, 1000), c(0, 500))) }
  if(k == 7) { ppp.dat <- ppp(x = raw_census7$PX, y = raw_census7$PY, marks = raw_census7$Latin, window = owin(c(0, 1000), c(0, 500))) }
  if(k == 8) { ppp.dat <- ppp(x = raw_census8$PX, y = raw_census8$PY, marks = raw_census8$Latin, window = owin(c(0, 1000), c(0, 500))) }
  
  bci_spatial_bc <- data.frame(Year = census_years[k], Quadrat_size = quadrat_sizes, Mean_quadrat_dissimilarity = NA, SD_quadrat_dissimilarity = NA)
  
  # Perform analyses for each quadrat size and census combination
  for(i in 1:length(quadrat_sizes)) { # quadrat size loop
    quad.df <- data.frame(x = ppp.dat$x, y = ppp.dat$y, species = ppp.dat$marks, year = census_years[k])
    quad.df$xt <- cut(quad.df$x, seq(xmin, xmax, quadrat_sizes[i]))
    quad.df$yt <- cut(quad.df$y, seq(ymin, ymax, quadrat_sizes[i]))
    quad.df <- quad.df %>% 
      separate(xt, into = c("temp", "temp1"), sep = ",", remove = F) %>% 
      separate(temp, 1, into = c("temp2", "qx")) %>%
      separate(yt, into = c("temp3", "temp4"), sep = ",", remove = F) %>% 
      separate(temp3, 1, into = c("temp5", "qy")) %>%
      mutate(qx = as.numeric(qx), qy = as.numeric(qy)) %>%
      dplyr::select(year, qx, qy, x, y, species, xt, yt) 
      
    quad.df$qx <- ifelse(quad.df$x == 0, 0, quad.df$qx) # where stems occur on the boundary, correctly assign coords
    quad.df$qy <- ifelse(quad.df$y == 0, 0, quad.df$qy)
    quad.df$Quadrat_size <- quadrat_sizes[i]
    
      quad.df <- quad.df %>%
        unite(quadrat, c(qx, qy), sep = "_", remove = F) %>%
        unite(year_quadrat, c(year, qx, qy), sep = "_", remove = F) %>%
        arrange(year_quadrat, species)
    
    # Generate datasets in wide format
    wide_quad.df <- quad.df %>% 
      mutate(Presence = 1) %>%  # add a column to indicate presence
      group_by(quadrat, species) %>% # for each species at each census and quadrat...
      summarise(Abundance = sum(Presence)) %>% # ...sum the number of individuals 
      spread(key = species, value = Abundance, fill = 0) %>% # create wide dataset using species names and presence
      ungroup() %>%
      dplyr::select(order(colnames(.)), -quadrat) %>% # sort dataset by species alphabetically 
      as.data.frame()
    rownames(wide_quad.df) <- unique(quad.df$quadrat) 
    
    # Calculate mean spatial quadrat dissimilarity (among quadrats within censuses)
    bci_spatial_bc$Mean_quadrat_dissimilarity[i] <- mean(vegdist(wide_quad.df, method = "bray"))
    bci_spatial_bc$SD_quadrat_dissimilarity[i] <- sd(vegdist(wide_quad.df, method = "bray"))
    
    # Rename datasets for output
    quad.df$presence <- 1
    nam <- paste("bci_census", k, "quad", quadrat_sizes[i], "long", sep = "_")
    assign(nam, quad.df)
     
    # nam <- paste("bci_census", k, "quad", quadrat_sizes[i], "wide", sep = "_")
    # assign(nam, wide_quad.df)
    
    if(i == 1) { bci_quad_long <- quad.df }
    if(i > 0)  { bci_quad_long <- rbind(bci_quad_long, quad.df) }
    
  } # end i loop (quadrat sizes)
  
  if(k == 1) { bci_spatial_bc_result <- bci_spatial_bc 
               bci_quad_long_out <- bci_quad_long }
  if(k > 1) { bci_spatial_bc_result <- rbind(bci_spatial_bc_result, bci_spatial_bc) 
              bci_quad_long_out <- rbind(bci_quad_long_out, bci_quad_long) }

} # end k loop (censuses)

#######################################################################################
# Use raw dissimilarities in space and time to determine if quadrats are diverging over time
#######################################################################################

# Calculate mean dissimilarity among quadrats at each census time (for quadrats of different sizes)
bci_spatial_bc_result %>% 
  mutate(Census = as.factor(Year)) %>%
  ggplot(aes(x = Quadrat_size, y = Mean_quadrat_dissimilarity, colour = Census)) +
    geom_line(size = 1, alpha = 0.8) +
  geom_errorbar(aes(ymin = Mean_quadrat_dissimilarity - SD_quadrat_dissimilarity, ymax = Mean_quadrat_dissimilarity + SD_quadrat_dissimilarity), position = "dodge", width = 0.25) +
    xlab("Quadrat size (m)") + ylab("Mean quadrat Bray-Curtis dissimlarity\n(SD error bars)") +
    theme_bw(base_size = 14)
ggsave("BCI_temporal_dissimilarity_quad_size.png", height = 6, width = 12, dpi = 600)

# Calculate mean temporal quadrat dissimilarity (among censuses for quadrats)
for(i in 1: length(quadrat_sizes)){ # loop through each quadrat size
  temp.df <- bci_quad_long_out[bci_quad_long_out$Quadrat_size == quadrat_sizes[i], ] # extract data for each quad size
  quadrat.list <- unique(temp.df$quadrat)
  
  bci_temporal_bc <- data.frame(Quadrat_size = quadrat_sizes[i], Quadrat = quadrat.list, Mean_census_dissimilarity = NA, SD_census_dissimilarity = NA) # create output object
  
  for(q in 1: length(unique(temp.df$quadrat))) { # loop through each quadrat
    wide_quad.df <- temp.df[temp.df$quadrat == quadrat.list[q], ] %>% # extract data for each quadrat and make wide format df
      mutate(Presence = 1) %>%  # add a column to indicate presence
      group_by(year, species) %>% # for each species at each census and quadrat...
      summarise(Abundance = sum(Presence)) %>% # ...sum the number of individuals 
      spread(key = species, value = Abundance, fill = 0) %>% # create wide dataset using species names and presence
      ungroup() %>%
      dplyr::select(order(colnames(.)), -year) %>% # sort dataset by species alphabetically 
      as.data.frame() 
    
    # Calculate mean tempoarl quadrat dissimilarity (among censuses for each quadrat)
    bci_temporal_bc$Mean_census_dissimilarity[q] <- mean(vegdist(wide_quad.df, method = "bray"))
    bci_temporal_bc$SD_census_dissimilarity[q] <- sd(vegdist(wide_quad.df, method = "bray"))
    
    } # end q loop through each quadrat

   # add results to output object
   if(i == 1) { bci_temporal_bc_result <- bci_temporal_bc }
   if(i > 1) { bci_temporal_bc_result <- rbind(bci_temporal_bc_result, bci_temporal_bc) }
  
 } # end i loop through each quadrat size

# Plot heat maps of temporal dissimilarity by quadrats for a range of quadrat sizes 
bci_temporal_bc_result %>%
  separate(Quadrat, into = c("x", "y"), sep = "_") %>%
  mutate(x = (as.numeric(x) + Quadrat_size/2), y = (as.numeric(y) + Quadrat_size/2)) %>%
  ggplot(aes(x = x, y = y, fill = Mean_census_dissimilarity, colour = Mean_census_dissimilarity)) +
  geom_tile(aes(fill = Mean_census_dissimilarity, width = Quadrat_size, height = Quadrat_size), colour = "grey50", hjust = 0, vjust = 0) +
  scale_fill_gradientn(colours = terrain.colors(10)) +
  coord_fixed(ratio = 1) +
  xlim(c(0, 1000)) + ylim(c(0, 500)) +
  xlab("X (m)") + ylab("Y (m)") +
  facet_grid(Quadrat_size ~ .) +
  theme_bw(base_size = 14)
ggsave("BCI_temporal_dissimilarity_quad_heatmap.png", height = 6, width = 6, dpi = 600)

bci_temporal_bc_result %>%
  separate(Quadrat, into = c("x", "y"), sep = "_") %>%
  mutate(x = (as.numeric(x) + Quadrat_size/2), y = (as.numeric(y) + Quadrat_size/2)) %>%
  ggplot(aes(x = x, y = y, fill = SD_census_dissimilarity, colour = SD_census_dissimilarity)) +
  geom_tile(aes(fill = SD_census_dissimilarity, width = Quadrat_size, height = Quadrat_size), colour = "grey50", hjust = 0, vjust = 0) +
  coord_fixed(ratio = 1) +
  xlim(c(0, 1000)) + ylim(c(0, 500)) +
  xlab("X (m)") + ylab("Y (m)") +
  facet_grid(Quadrat_size ~ .) +
  theme_bw(base_size = 14)
ggsave("BCI_temporal_dissimilarity_SD_quad_heatmap.png", height = 6, width = 6, dpi = 600)


##########################################################
# Zeta diversity
##########################################################

bci_20$quadrat_size <- 20
bci_50$quadrat_size <- 50
bci_125$quadrat_size <- 125
bci_250$quadrat_size <- 250
bci_500$quadrat_size <- 500

bci_quad_long <- rbind(bci_20, bci_50, bci_125, bci_250, bci_500)

quadrat_sizes <- c(20, 50, 125, 250, 500)

for(i in 1:length(unique(bci_quad_long$quadrat_size))) { # loop through the different quadrat sizes
  temp.df <- bci_quad_long[bci_quad_long$quadrat_size == quadrat_sizes[i], ] # subset each quadrat size
  quadrat.list <- unique(temp.df$quadrat)
  temp.df$presence <- 1
  
  for(q in 1: length(quadrat.list)) { # loop through each quadrat
    wide_quad.df <- temp.df[temp.df$quadrat == quadrat.list[q], ] %>% # extract data for each quadrat
      dplyr::select(year, species, presence) %>% # select required variables
      spread(key = species, value = presence, fill = 0) %>% # create presence absence matrix
      as.data.frame() # convert to a data frame object

    bci_zeta <- data.frame(Quadrat = NA, Quadrat_size = NA, Order = 1:7, Zeta_diversity = NA, Zeta_diversity_SD = NA)
    
    for(z in 1:7) { # loop through the orders of zeta
      bci_zeta$Quadrat[z] <- quadrat.list[q]
      bci_zeta$Quadrat_size[z] <- quadrat_sizes[i]
      
      # Calculate temporal normalised mean zeta diversity for each quadrat: Number of shared species among different numbers of samples (orders of zeta) divided by the minimum number of species in the sites of each specific sample 
      zeta.out <- Zeta.order.mc(wide_quad.df[, -1], order = z, normalize = "Simpson")
      bci_zeta$Order[z] <- zeta.out$zeta.order
      bci_zeta$Zeta_diversity[z] <- zeta.out$zeta.val
      bci_zeta$Zeta_diversity_SD[z] <- zeta.out$zeta.val.sd
      
    } # end z loop through each order of zeta

    if(q == 1) { bci_zeta_quad <- bci_zeta }
    if(q > 1) { bci_zeta_quad <- rbind(bci_zeta_quad, bci_zeta) }
    
  } # end q loop through each quadrat
  
  if(i == 1) { bci_zeta_result <- bci_zeta_quad }
  if(i > 1) { bci_zeta_result <- rbind(bci_zeta_result, bci_zeta_quad) }
  
} # end i loop through each quadrat size

### Plot zeta diversity results
bci_zeta_result %>%
  filter(Order != 1) %>%
  separate(Quadrat, into = c("x", "y"), sep = "_") %>%
  mutate(x = (as.numeric(x) + Quadrat_size/2), y = (as.numeric(y) + Quadrat_size/2)) %>%
  mutate(inv_Zeta = 1 - Zeta_diversity) %>%
  ggplot(aes(x = x, y = y, fill = inv_Zeta, colour = inv_Zeta)) +
  geom_tile(aes(fill = inv_Zeta, width = Quadrat_size, height = Quadrat_size), colour = "grey50", hjust = 0, vjust = 0) +
  scale_fill_gradientn(colours = terrain.colors(10)) +
  coord_fixed(ratio = 1) +
  xlim(c(0, 1000)) + ylim(c(0, 500)) +
  xlab("X (m)") + ylab("Y (m)") +
  facet_grid(Quadrat_size ~ Order) +
  theme_bw(base_size = 14)
ggsave("BCI_zeta_quad_heatmap.png", height = 8, width = 16, dpi = 600)

### Calculate observed inverse zeta values and their significance compared to random expectation 
for(i in 1:length(unique(bci_quad_long$quadrat_size))) { # loop through the different quadrat sizes
  temp.df <- bci_quad_long[bci_quad_long$quadrat_size == quadrat_sizes[i], ] # subset each quadrat size
  quadrat.list <- unique(temp.df$quadrat)
  temp.df$presence <- 1
  
  # Create output object for each quadrat
  bci_zeta <- data.frame(Quadrat = rep(NA, length(quadrat.list)), Quadrat_size = NA, Inv_zeta_2 = NA, Inv_zeta_3 = NA, Inv_zeta_4 = NA, Inv_zeta_5 = NA, Inv_zeta_6 = NA, Inv_zeta_7 = NA, mean_pred_inv_zeta_2 = NA, mean_pred_inv_zeta_3 = NA, mean_pred_inv_zeta_4 = NA, mean_pred_inv_zeta_5 = NA, mean_pred_inv_zeta_6 = NA, mean_pred_inv_zeta_7 = NA, p_zeta_2 = NA, p_zeta_3 = NA, p_zeta_4 = NA, p_zeta_5 = NA, p_zeta_6 = NA, p_zeta_7 = NA)
  
  for(q in 1: length(quadrat.list)) { # loop through each quadrat
    wide_quad.df <- temp.df[temp.df$quadrat == quadrat.list[q], ] %>% # extract data for each quadrat
      dplyr::select(year, species, presence) %>% # select required variables
      spread(key = species, value = presence, fill = 0) %>% # create presence absence matrix
      as.data.frame() # convert to a data frame object
    
      bci_zeta$Quadrat[q] <- quadrat.list[q] # add quadrat information
      bci_zeta$Quadrat_size[q] <- quadrat_sizes[i]
      
      # Calculate temporal normalised mean zeta diversity for each quadrat: Number of shared species among different numbers of samples (orders of zeta) divided by the minimum number of species in the sites of each specific sample 
      bci_zeta$Inv_zeta_2[q] <- 1 - Zeta.order.mc(wide_quad.df[, -1], order = 2, normalize = "Simpson")$zeta.val
      bci_zeta$Inv_zeta_3[q] <- 1 - Zeta.order.mc(wide_quad.df[, -1], order = 3, normalize = "Simpson")$zeta.val
      bci_zeta$Inv_zeta_4[q] <- 1 - Zeta.order.mc(wide_quad.df[, -1], order = 4, normalize = "Simpson")$zeta.val
      bci_zeta$Inv_zeta_5[q] <- 1 - Zeta.order.mc(wide_quad.df[, -1], order = 5, normalize = "Simpson")$zeta.val
      bci_zeta$Inv_zeta_6[q] <- 1 - Zeta.order.mc(wide_quad.df[, -1], order = 6, normalize = "Simpson")$zeta.val
      bci_zeta$Inv_zeta_7[q] <- 1 - Zeta.order.mc(wide_quad.df[, -1], order = 7, normalize = "Simpson")$zeta.val
      
      # Permute the community matrix in time only (species relative commonness and rarity is retained; only their presence at the different census times is permuted)
      bci.perm <- permatfull(wide_quad.df[-1], fixedmar = "columns", times = 199)
      
      # Create output object for permutations
      bci.perm.results <- data.frame(permutation = 1:length(bci.perm$perm), pred_inv_zeta_2 = NA, pred_inv_zeta_3 = NA, pred_inv_zeta_4 = NA, pred_inv_zeta_5 = NA, pred_inv_zeta_6 = NA, pred_inv_zeta_7 = NA)
      
      # loop through each permutation
      for (p in 1:length(bci.perm.results$pred_inv_zeta_2)){
        perm.df <- as.data.frame(bci.perm$perm[p][[1]])
        bci.perm.results$pred_inv_zeta_2[p] <- 1 - Zeta.order.mc(perm.df, order = 2, normalize = "Simpson")$zeta.val
        bci.perm.results$pred_inv_zeta_3[p] <- 1 - Zeta.order.mc(perm.df, order = 3, normalize = "Simpson")$zeta.val
        bci.perm.results$pred_inv_zeta_4[p] <- 1 - Zeta.order.mc(perm.df, order = 4, normalize = "Simpson")$zeta.val
        bci.perm.results$pred_inv_zeta_5[p] <- 1 - Zeta.order.mc(perm.df, order = 5, normalize = "Simpson")$zeta.val
        bci.perm.results$pred_inv_zeta_6[p] <- 1 - Zeta.order.mc(perm.df, order = 6, normalize = "Simpson")$zeta.val
        bci.perm.results$pred_inv_zeta_7[p] <- 1 - Zeta.order.mc(perm.df, order = 7, normalize = "Simpson")$zeta.val
        
        print(paste("Quadrat Size =", quadrat_sizes[i], "Quadrat =", quadrat.list[q], "Permutation =", p))
        
        } # end permutation loop
      
      # Store mean inverse zeta of the permuted values
      bci_zeta$mean_pred_inv_zeta_2[q] <- mean(bci.perm.results$pred_inv_zeta_2)
      bci_zeta$mean_pred_inv_zeta_3[q] <- mean(bci.perm.results$pred_inv_zeta_3)
      bci_zeta$mean_pred_inv_zeta_4[q] <- mean(bci.perm.results$pred_inv_zeta_4)
      bci_zeta$mean_pred_inv_zeta_5[q] <- mean(bci.perm.results$pred_inv_zeta_5)
      bci_zeta$mean_pred_inv_zeta_6[q] <- mean(bci.perm.results$pred_inv_zeta_6)
      bci_zeta$mean_pred_inv_zeta_7[q] <- mean(bci.perm.results$pred_inv_zeta_7)
      
      # Generate p-values for comparing observed to predicted inverse zeta
      bci_zeta$p_zeta_2[q] <- min((sum(bci.perm.results$pred_inv_zeta_2 < bci_zeta$Inv_zeta_2[q]) + 1), (sum(bci.perm.results$pred_inv_zeta_2 >= bci_zeta$Inv_zeta_2[q]) + 1))*2 / (length(bci.perm.results$pred_inv_zeta_2) + 1)
      bci_zeta$p_zeta_3[q] <- min((sum(bci.perm.results$pred_inv_zeta_3 < bci_zeta$Inv_zeta_3[q]) + 1), (sum(bci.perm.results$pred_inv_zeta_3 >= bci_zeta$Inv_zeta_3[q]) + 1))*2 / (length(bci.perm.results$pred_inv_zeta_3) + 1)
      bci_zeta$p_zeta_4[q] <- min((sum(bci.perm.results$pred_inv_zeta_4 < bci_zeta$Inv_zeta_4[q]) + 1), (sum(bci.perm.results$pred_inv_zeta_4 >= bci_zeta$Inv_zeta_4[q]) + 1))*2 / (length(bci.perm.results$pred_inv_zeta_4) + 1)
      bci_zeta$p_zeta_5[q] <- min((sum(bci.perm.results$pred_inv_zeta_5 < bci_zeta$Inv_zeta_5[q]) + 1), (sum(bci.perm.results$pred_inv_zeta_5 >= bci_zeta$Inv_zeta_5[q]) + 1))*2 / (length(bci.perm.results$pred_inv_zeta_5) + 1)
      bci_zeta$p_zeta_6[q] <- min((sum(bci.perm.results$pred_inv_zeta_6 < bci_zeta$Inv_zeta_6[q]) + 1), (sum(bci.perm.results$pred_inv_zeta_6 >= bci_zeta$Inv_zeta_6[q]) + 1))*2 / (length(bci.perm.results$pred_inv_zeta_6) + 1)
      bci_zeta$p_zeta_7[q] <- min((sum(bci.perm.results$pred_inv_zeta_7 < bci_zeta$Inv_zeta_7[q]) + 1), (sum(bci.perm.results$pred_inv_zeta_7 >= bci_zeta$Inv_zeta_7[q]) + 1))*2 / (length(bci.perm.results$pred_inv_zeta_7) + 1)
      
      #head(bci_zeta)
      
  } # end q loop through each quadrat
  
  if(i == 1) { bci_zeta_result <- bci_zeta }
  if(i > 1) { bci_zeta_result <- rbind(bci_zeta_result, bci_zeta) }
  
} # end i loop through each quadrat size

### Tidy dataset
bci_obs_zeta <- bci_zeta_result %>%
  dplyr::select(Quadrat_size, Quadrat, Inv_zeta_2, Inv_zeta_3, Inv_zeta_4, Inv_zeta_5, Inv_zeta_6, Inv_zeta_7) %>%
  gather(key = order, value = Inv_zeta, -Quadrat_size, -Quadrat) %>%
  separate(order, sep = -1, into = c("temp1", "Order")) %>%
  dplyr::select(Quadrat_size, Quadrat, Order, Inv_zeta) %>% 
  separate(Quadrat, into = c("x", "y"), sep = "_") %>%
  mutate(x = (as.numeric(x) + Quadrat_size/2), y = (as.numeric(y) + Quadrat_size/2))

bci_pred_zeta <- bci_zeta_result %>%
  dplyr::select(Quadrat_size, Quadrat, mean_pred_inv_zeta_2, mean_pred_inv_zeta_3, mean_pred_inv_zeta_4, mean_pred_inv_zeta_5, mean_pred_inv_zeta_6, mean_pred_inv_zeta_7) %>%
  gather(key = order, value = Mean_pred_inv_zeta, -Quadrat_size, -Quadrat) %>%
  dplyr::select(Mean_pred_inv_zeta)

bci_pvalue <- bci_zeta_result %>%
  dplyr::select(Quadrat_size, Quadrat, p_zeta_2, p_zeta_3, p_zeta_4, p_zeta_5, p_zeta_6, p_zeta_7) %>%
  gather(key = order, value = P_value, -Quadrat_size, -Quadrat) %>%
  dplyr::select(P_value)
  
bci_zeta_gdat <- cbind(bci_obs_zeta, bci_pred_zeta, bci_pvalue)

### Plot zeta diversity results

# Observed
bci_zeta_gdat %>%
  ggplot(aes(x = x, y = y, fill = Inv_zeta, colour = Inv_zeta)) +
  geom_tile(aes(fill = Inv_zeta, width = Quadrat_size, height = Quadrat_size), colour = "grey50", hjust = 0, vjust = 0) +
  scale_fill_gradientn(colours = terrain.colors(10), name = "1 - Zeta") +
  coord_fixed(ratio = 1) +
  xlim(c(0, 1000)) + ylim(c(0, 500)) +
  xlab("X (m)") + ylab("Y (m)") +
  facet_grid(Quadrat_size ~ Order) +
  theme_bw(base_size = 18)
ggsave("BCI_zeta_quad_heatmap.png", height = 8, width = 16, dpi = 600)

# Observed-Expected
bci_zeta_gdat %>%
  mutate(OE_zeta = Inv_zeta - Mean_pred_inv_zeta) %>%
  ggplot(aes(x = x, y = y, fill = OE_zeta, colour = OE_zeta)) +
  geom_tile(aes(fill = OE_zeta, width = Quadrat_size, height = Quadrat_size), colour = "grey50", hjust = 0, vjust = 0) +
  scale_fill_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0, space = "Lab", na.value = "grey50", guide = "colourbar", name = "Obs. - Exp.") +
  coord_fixed(ratio = 1) +
  xlim(c(0, 1000)) + ylim(c(0, 500)) +
  xlab("X (m)") + ylab("Y (m)") +
  facet_grid(Quadrat_size ~ Order) +
  theme_bw(base_size = 18)
ggsave("BCI_OE_zeta_quad_heatmap.png", height = 8, width = 16, dpi = 600)

# P-value
bci_zeta_gdat %>%
  ggplot(aes(x = x, y = y, fill = P_value, colour = P_value)) +
  geom_tile(aes(fill = P_value, width = Quadrat_size, height = Quadrat_size), colour = "grey50", hjust = 0, vjust = 0) +
  scale_fill_gradient(low = "red", high = "white", space = "Lab", na.value = "grey50", guide = "colourbar", limits=c(0,1), name = "P value") +
  coord_fixed(ratio = 1) +
  xlim(c(0, 1000)) + ylim(c(0, 500)) +
  xlab("X (m)") + ylab("Y (m)") +
  facet_grid(Quadrat_size ~ Order) +
  theme_bw(base_size = 18)
ggsave("BCI_P_zeta_quad_heatmap.png", height = 8, width = 16, dpi = 600)

# Significance
bci_zeta_gdat %>%
  mutate(Significance = factor(ifelse(P_value < 0.05,"Sig.","Non-sig."))) %>%
  ggplot(aes(x = x, y = y, fill = Significance, colour = Significance)) +
  geom_tile(aes(fill = Significance, width = Quadrat_size, height = Quadrat_size), colour = "grey50", hjust = 0, vjust = 0) +
  scale_fill_manual(values = c("steelblue3", "firebrick3")) +
  coord_fixed(ratio = 1) +
  xlim(c(0, 1000)) + ylim(c(0, 500)) +
  xlab("X (m)") + ylab("Y (m)") +
  facet_grid(Quadrat_size ~ Order) +
  theme_bw(base_size = 18)
ggsave("BCI_Sig_quad_heatmap.png", height = 8, width = 16, dpi = 600)

##########################################################
# Turnover
##########################################################

# One of the commonly used methods is that of Diamond (1969) where the percentage turnover rate, T, is calculated as T = 100 ?? (E + C) / (S1 + S2), where E is the number of species that went extinct between two time points, C is the number of species that colonised between the same two time points, S1 is the number of species present at time 1 and S2 is the number of species present at time 2

# Here we use the 'turnover' function in the codyn package which calculates turnover as the proportion of species either gained or lost relative to the total number of species observed across both time periods.

for(i in 1:length(unique(bci_quad_long$quadrat_size))) { # loop through the different quadrat sizes
  temp.df <- bci_quad_long[bci_quad_long$quadrat_size == quadrat_sizes[i], ] %>% # subset quadrat sizes
    dplyr::select(year, species, quadrat, abundance) 

  bci_obs_turnover_out <- turnover(temp.df, time.var = "year", species.var = "species", abundance.var = "abundance", replicate.var = "quadrat", metric = "total") %>%
    mutate(Quadrat_size = quadrat_sizes[i]) %>% # add quadrat size 
    rename(Obs_Turnover = total) 
  
  # Permute the community matrix in time only (species relative commonness and rarity is retained; only their presence at the different census times is permuted)
  
  quadrat.list <- unique(temp.df$quadrat) # make a list of quadrat names
  
  for(q in 1: length(quadrat.list)) { # loop through each quadrat
    wide_quad.df <- temp.df[temp.df$quadrat == quadrat.list[q], ] %>% # extract data for each quadrat
      dplyr::select(year, species, abundance) %>% # select required variables
      spread(key = species, value = abundance, fill = 0) %>% # create abundance matrix
      as.data.frame() # convert to a data frame object
    
  bci.perm <- permatfull(wide_quad.df[-1], fixedmar = "columns", times = 199)
  
  # loop through each permutation
  for (p in 1:length(bci.perm$perm)) { # loop through the permutations
    perm.df <- as.data.frame(bci.perm$perm[[p]]) %>%
      mutate(year = c(1982, 1985, 1990, 1995, 2000, 2005, 2010, 2015)) %>%
      gather(species, abundance, -year)
    
    bci_pred_turnover_quad <- as.data.frame(turnover(perm.df, time.var = "year", species.var = "species", abundance.var = "abundance", metric = "total"))
    
    bci_pred_turnover_quad$Permutation <- p # add permutation number
    bci_pred_turnover_quad$Quadrat <- quadrat.list[q] # add quadrat name
    bci_pred_turnover_quad$Quadrat_size <- quadrat_sizes[i] # add quadrat size
    
    print(paste("Quadrat Size =", quadrat_sizes[i], "Quadrat =", quadrat.list[q], "Permutation =", p))
    
    if(p == 1) { out1 <- bci_pred_turnover_quad }
    if(p > 1) { out1 <- rbind(out1, bci_pred_turnover_quad) }
    
  } # end p loop through permutations
  
  if(q == 1) { out2 <- out1 }
  if(q > 1) { out2 <- rbind(out2, out1) }
  
  } # end q loop through quadrat list
  
  if(i == 1) { bci_turnover_perm_out <- out2
               bci_turnover_obs_out <- bci_obs_turnover_out }
  if(i > 1) {  bci_turnover_perm_out <- rbind(bci_turnover_perm_out, out2) 
               bci_turnover_obs_out <- rbind(bci_turnover_obs_out, bci_obs_turnover_out)  }
  
  } # end i loop through quadrat sizes
  
head(bci_turnover_perm_out)
head(bci_turnover_obs_out)

# Compare predicted distributions to the observed turnover values for each quadrat and time period

years <- unique(bci_turnover_perm_out$year)

for(i in 1:length(unique(bci_quad_long$quadrat_size))) { # loop through the different quadrat sizes
  temp.df <- bci_turnover_perm_out[bci_turnover_perm_out$Quadrat_size == quadrat_sizes[i], ] 
  quadrat.list <- unique(temp.df$Quadrat) # make a list of quadrat names
  
  for(q in 1: length(quadrat.list)) { # loop through each quadrat
    temp.quad.df <- temp.df[temp.df$Quadrat == quadrat.list[q], ] 
    
    for(t in 1:length(unique(temp.quad.df$year))){
      temp.year.df <- temp.quad.df[temp.quad.df$year == years[t], ]
      obs.turnover <- bci_turnover_obs_out %>% dplyr::filter(year == years[t], quadrat == quadrat.list[q], Quadrat_size == quadrat_sizes[i]) %>% dplyr::select(Obs_Turnover) %>% as.numeric(.) # obtain correct observed value
      p_value <-  min((sum(temp.year.df$total < obs.turnover) + 1), (sum(temp.year.df$total >= obs.turnover) + 1))*2 / (length(temp.year.df$total) + 1) # calculate p-value
      mean_pred_turnover <- mean(temp.year.df$total) # calculate mean predicted value
      o_e_turnover <- obs.turnover - mean_pred_turnover # calculate observed minus expected value
      significance <- ifelse(p_value > 0.05, "Non-sig.", "Sig.") # categorise significance
      temp_out <- data.frame(Quadrat_size = quadrat_sizes[i], year = years[t], Quadrat = quadrat.list[q], Obs_turnover = obs.turnover, Mean_pred_turnover = mean_pred_turnover, Obs_exp_turnover = o_e_turnover, P_value = p_value, Significance = significance) # write to output dataframe
      
      if(t == 1) { out1 <- temp_out }
      if(t > 1) { out1 <- rbind(out1, temp_out) }
      
      print(paste("Quadrat Size =", quadrat_sizes[i], "Quadrat =", quadrat.list[q], "Year =", years[t]))
      
  } # end t loop through years
    
    if(q == 1) { out2 <- out1 }
    if(q > 1) { out2 <- rbind(out2, out1) }
    
  } # end q loop through quadrat list

  if(i == 1) { bci_turnover_results <- out2 }
  if(i > 1) { bci_turnover_results <- rbind(bci_turnover_results, out2) }
  
    } # end i loop through quadrat sizes

bci_gdat <- bci_turnover_results %>%
  separate(Quadrat, into = c("x", "y"), sep = "_") %>%
  mutate(x = (as.numeric(x) + Quadrat_size/2), y = (as.numeric(y) + Quadrat_size/2)) %>%
  mutate(year = as.factor(year)) %>%
  mutate(Time_period = fct_recode(year, "1982 - 1985" = "1985", "1985 - 1990" = "1990", "1990 - 1995" = "1995", "1995 - 2000" = "2000", "2000 - 2005" = "2005", "2005 - 2010" = "2010", "2010 - 2015" = "2015"))
  
### Plot turnover results

# Observed
bci_gdat %>%
  ggplot(aes(x = x, y = y, fill = Obs_turnover, colour = Obs_turnover)) +
  geom_tile(aes(fill = Obs_turnover, width = Quadrat_size, height = Quadrat_size), colour = "grey50", hjust = 0, vjust = 0) +
  scale_fill_gradientn(colours = terrain.colors(10), name = "Turnover") +
  coord_fixed(ratio = 1) +
  xlim(c(0, 1000)) + ylim(c(0, 500)) +
  xlab("X (m)") + ylab("Y (m)") +
  facet_grid(Quadrat_size ~ Time_period) +
  theme_bw(base_size = 18)
ggsave("BCI_obs_turnover_quad_heatmap.png", height = 8, width = 16, dpi = 600)

# Observed-Expected
bci_gdat %>%
  ggplot(aes(x = x, y = y, fill = Obs_exp_turnover, colour = Obs_exp_turnover)) +
  geom_tile(aes(fill = Obs_exp_turnover, width = Quadrat_size, height = Quadrat_size), colour = "grey50", hjust = 0, vjust = 0) +
  scale_fill_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0, space = "Lab", na.value = "grey50", guide = "colourbar", name = "Obs. - Exp.") +
  coord_fixed(ratio = 1) +
  xlim(c(0, 1000)) + ylim(c(0, 500)) +
  xlab("X (m)") + ylab("Y (m)") +
  facet_grid(Quadrat_size ~ Time_period) +
  theme_bw(base_size = 18)
ggsave("BCI_OE_turnover_quad_heatmap.png", height = 8, width = 16, dpi = 600)

# P-value
bci_gdat %>%
  ggplot(aes(x = x, y = y, fill = P_value, colour = P_value)) +
  geom_tile(aes(fill = P_value, width = Quadrat_size, height = Quadrat_size), colour = "grey50", hjust = 0, vjust = 0) +
  scale_fill_gradient(low = "red", high = "white", space = "Lab", na.value = "grey50", guide = "colourbar", limits=c(0,1), name = "P value") +
  coord_fixed(ratio = 1) +
  xlim(c(0, 1000)) + ylim(c(0, 500)) +
  xlab("X (m)") + ylab("Y (m)") +
  facet_grid(Quadrat_size ~ Time_period) +
  theme_bw(base_size = 18)
ggsave("BCI_P_turnover_quad_heatmap.png", height = 8, width = 16, dpi = 600)

# Significance
bci_gdat %>%
  mutate(Significance = factor(ifelse(P_value < 0.05,"Sig.","Non-sig."))) %>%
  ggplot(aes(x = x, y = y, fill = Significance, colour = Significance)) +
  geom_tile(aes(fill = Significance, width = Quadrat_size, height = Quadrat_size), colour = "grey50", hjust = 0, vjust = 0) +
  scale_fill_manual(values = c("steelblue3", "firebrick3")) +
  coord_fixed(ratio = 1) +
  xlim(c(0, 1000)) + ylim(c(0, 500)) +
  xlab("X (m)") + ylab("Y (m)") +
  facet_grid(Quadrat_size ~ Time_period) +
  theme_bw(base_size = 18)
ggsave("BCI_Sig_turnover_heatmap.png", height = 8, width = 16, dpi = 600)

#save.image("BCI_output_18-05-13.RData")
#load("BCI_output_18-05-13.RData")