## Fig_S3 in supporting information

#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)
})

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

#Bioclimatic-zone specific temporal responses
species_temp = c('Moth')
data_temp = data.frame(Species = species_temp,
                       Response = factor(rep(c('CTI', 'CTIdiv' ,'CTBI', 'CTBIdiv', 'Warm', 'Cold'), each=5), levels = c('CTI', 'CTIdiv' ,'CTBI', 'CTBIdiv', 'Cold', 'Warm')),
                       Area = factor(rep(c('All sites', 'Hemiboreal', 'Southern boreal', 'Middle boreal', 'Northern boreal'), 6),
                                     levels = c('All sites', 'Hemiboreal', 'Southern boreal', 'Middle boreal', 'Northern boreal')),
                       Mean = NA, Min = NA, Max = NA)

#CTI, thermal range and thermal diversity
# Mean
resp_temp = c(models[[1]]$summary.lincomb.derived[1:5,2],
              models[[2]]$summary.lincomb.derived[1:5,2],
              models[[3]]$summary.lincomb.derived[1:5,2],
              models[[4]]$summary.lincomb.derived[1:5,2],
              models[[5]]$summary.lincomb.derived[1:5,2],
              models[[6]]$summary.lincomb.derived[1:5,2])

data_temp$Mean = resp_temp

# Min
resp_temp = c(models[[1]]$summary.lincomb.derived[1:5,4],
              models[[2]]$summary.lincomb.derived[1:5,4],
              models[[3]]$summary.lincomb.derived[1:5,4],
              models[[4]]$summary.lincomb.derived[1:5,4],
              models[[5]]$summary.lincomb.derived[1:5,4],
              models[[6]]$summary.lincomb.derived[1:5,4])
data_temp$Min = resp_temp

# Max
resp_temp = c(models[[1]]$summary.lincomb.derived[1:5,6],
              models[[2]]$summary.lincomb.derived[1:5,6],
              models[[3]]$summary.lincomb.derived[1:5,6],
              models[[4]]$summary.lincomb.derived[1:5,6],
              models[[5]]$summary.lincomb.derived[1:5,6],
              models[[6]]$summary.lincomb.derived[1:5,6])
data_temp$Max = resp_temp


#panel letters
graphLabels = data.frame(Response = factor(c("CTI", "CTIdiv", "CTBI", "CTBIdiv"), levels = c("CTI", "CTIdiv", "CTBI", "CTBIdiv")),
                         Species = .5, Mean = c(0.0275, 0.0009, 0.002, 0.0001),
                          Letter = c("A)", "B)", "C)", "D)"))

g1 = ggplot(data = data_temp[data_temp$Response %in% c('CTI', 'CTIdiv' ,'CTBI', 'CTBIdiv'),], aes(x = Species, y = Mean, fill = Area)) +
  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) +
  geom_text(data = graphLabels, aes(x = Species, y = Mean, label = Letter), inherit.aes = FALSE) +
  theme_classic() +
  geom_hline(yintercept = 0) +
  facet_wrap(~Response, scales = "free") +
  ylab(expression(paste(Delta, ' °C / year'))) +
  xlab('') +
  scale_fill_moma_d("VanGogh", name = '') +
  theme(legend.position = 'bottom', legend.text=element_text(size=6), axis.text.x = element_blank()) +
  guides(fill = guide_legend(nrow = 2))


g2 = ggplot(data = data_temp[data_temp$Area=='All sites' & data_temp$Response %in% c('Cold', 'Warm'),], aes(x = Species, y = Mean, fill = Response)) +
  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) +
  geom_text(x = .5, y = .0125, label = 'E)') +
  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(.75,.2), axis.text.x = element_blank())

cm2in = 0.393700787

png(paste0('Figures/Fig_S3.png'), width = (210-2*2)*.8, height = (297-2.54*2)/2, res = 600, units = 'mm')
grid.arrange(g1,g2, nrow =1, widths = c(.6,.4))
dev.off()

# pdf(paste0('Figures/Fig2.pdf'), width = 11*cm2in, height = 7*cm2in)
# g_layout = rbind(c(1,2),c(1,3))
# grid.arrange(g1, g2, g3, layout_matrix = g_layout, widths = c(.8,.2))
# dev.off()
