##Figure 3 of the manuscript

#set path
try(setwd(''))

suppressPackageStartupMessages({
  library(terra)
  library(tidyverse)
  library(vegan)
  library(reshape)
  library(MASS)
  library(sf)
  library(rnaturalearth)
  library(rnaturalearthdata)
  library(paletteer)
  library(ggrepel)
  library(gridExtra)
  library(RColorBrewer)
  library(INLA)
  library(pracma)
  library(viridis)
  library(ggsci)
  library(wesanderson)
  library(betapart)
  library(parallel)
  library(stars)
  library(adegenet)
  library(spdep)
  library(MoMAColors)
  library(tidyterra)
})


#coordinate system EPSG: 3067
#get a background map
world = ne_countries(scale = "medium", returnclass = "sf")
world_plane = st_transform(world, 3067)
#bioclimatic zones
bioclim = st_read('data/bioclimatic_zones/BioclimaticZonesFinland.shp')

##read in model structures and data
load('data/Model_fits.RData')

#Bioclimatic-zone specific temporal responses
species_temp = c('Moth', 'Butterfly', 'Bird', 'Plant', 'Phytopl.')
area_temp = c('All sites', 'Hemiboreal', 'Southern boreal', 'Middle boreal', 'Northern boreal')

# species_list = c('moth', 'butterfly', 'bird', 'plant', 'phytopl.')
data_temp = data.frame(Species = factor(rep(species_temp, each = 5*2), levels = c('Plant', 'Phytopl.', 'Butterfly', 'Moth', 'Bird')),
                       STI = factor(rep(rep(c('Warm', 'Cold'), each = 5), 5), levels = c('Cold', 'Warm')),
                       Area = factor(rep(rep(c('All sites', 'Hemiboreal', 'Southern boreal', 'Middle boreal', 'Northern boreal'), 2), 5),
                                     levels = c('All sites', 'Hemiboreal', 'Southern boreal', 'Middle boreal', 'Northern boreal')),
                       Mean = NA, Min = NA, Max = NA)

# store t-test scores
t_test_score = {}
t_test_score_abs = {}

for (i_sp in 1:length(species_temp)) {
  print(i_sp)
  
  #CTI, thermal range and thermal diversity
  # Mean
  resp_temp = c(models[[i_sp]][[5]]$summary.lincomb.derived[1:5,2],
                models[[i_sp]][[6]]$summary.lincomb.derived[1:5,2])
  data_temp$Mean[data_temp$Species==species_temp[i_sp]] = resp_temp
  
  # Min
  resp_temp = c(models[[i_sp]][[5]]$summary.lincomb.derived[1:5,4],
                models[[i_sp]][[6]]$summary.lincomb.derived[1:5,4])
  data_temp$Min[data_temp$Species==species_temp[i_sp]] = resp_temp
  
  # Max
  resp_temp = c(models[[i_sp]][[5]]$summary.lincomb.derived[1:5,6],
                models[[i_sp]][[6]]$summary.lincomb.derived[1:5,6])
  data_temp$Max[data_temp$Species==species_temp[i_sp]] = resp_temp
  
  #t-test score for estimates
  mean_temp = c(models[[i_sp]][[5]]$summary.lincomb.derived[1,2],
                models[[i_sp]][[6]]$summary.lincomb.derived[1,2])
  sd_temp = c(models[[i_sp]][[5]]$summary.lincomb.derived[1,3],
              models[[i_sp]][[6]]$summary.lincomb.derived[1,3])
  sample_temp = nrow(xy_unique[[i_sp]])
  
  pooled_se = sqrt((sd_temp[1]^2 / sample_temp) + (sd_temp[2]^2 / sample_temp))
  t_stat = (mean_temp[1] - mean_temp[2]) / pooled_se
  
  df = ((sd_temp[1]^2 / sample_temp) + (sd_temp[2]^2 / sample_temp))^2 /
    (((sd_temp[1]^2 / sample_temp)^2 / (sample_temp - 1)) + ((sd_temp[2]^2 / sample_temp)^2 / (sample_temp - 1)))
  t_test_score[[i_sp]] = 2 * pt(-abs(t_stat), df)
  
  #t-test score for absolute estimates
  if (i_sp < 4) {
    mean_temp = c(models[[i_sp]][[5]]$summary.lincomb.derived[1,2],
                  -models[[i_sp]][[6]]$summary.lincomb.derived[1,2])
    sd_temp = c(models[[i_sp]][[5]]$summary.lincomb.derived[1,3],
                models[[i_sp]][[6]]$summary.lincomb.derived[1,3])
    sample_temp = nrow(xy_unique[[i_sp]])
    
    pooled_se = sqrt((sd_temp[1]^2 / sample_temp) + (sd_temp[2]^2 / sample_temp))
    t_stat = (mean_temp[1] - mean_temp[2]) / pooled_se
    
    df = ((sd_temp[1]^2 / sample_temp) + (sd_temp[2]^2 / sample_temp))^2 /
      (((sd_temp[1]^2 / sample_temp)^2 / (sample_temp - 1)) + ((sd_temp[2]^2 / sample_temp)^2 / (sample_temp - 1)))
    t_test_score_abs[[i_sp]] = 2 * pt(-abs(t_stat), df)
  }
  
}


# data_temp$Species = factor(data_temp$Species, levels = c('Plant', 'Phytoplankton', 'Butterfly', 'Moth', 'Bird'))
# data_temp$Area = factor(data_temp$Area)
# data_temp$STI = factor(data_temp$STI, levels = c('Cold', 'Cool', 'Warm', 'Hot'))

g1 = ggplot(data = data_temp[data_temp$Area=='All sites',], aes(x = Species, y = Mean, fill = STI)) +
  geom_bar(stat = 'identity', position=position_dodge(.8), width=.8) +
  geom_point(position=position_dodge(.8), size = .6) +
  geom_errorbar(aes(ymin = Min, ymax = Max), position=position_dodge(.8), width = .5) +
  theme_classic() +
  geom_hline(yintercept = 0) +
  ylab(expression(paste(Delta, ' log(abundance) / year'))) +
  xlab('') +
  scale_fill_manual(values = brewer.pal(n = 9, name = "RdYlBu")[c(7,3)], name = 'STI') +
  theme(legend.position = c(.2,.2))

cm2in = 0.393700787

png(paste0('Figures/Fig3.png'), width = (210-2*2)/1.5, height = (297-2.54*2)/3, res = 600, units = 'mm')
g1
dev.off()

tiff(paste0('Figures/Fig3.tif'), width = (210-2*2)/1.5, height = (297-2.54*2)/3, res = 200, units = 'mm')
g1
dev.off()

pdf(paste0('Figures/Fig3.pdf'), width = 11*cm2in, height = 9*cm2in)
g1
dev.off()

