code.R

# Code to accompany paper: Temperate functional niche availability not resident-invader competition shapes tropicalisation in reef fishes

# Paper authors: Mark Miller, James Reimer, Katie Cook, John M. Pandolfi, Brigitte Sommer, Masami Obuchi, Maria Beger
# Code author: Mark Miller mark.gr.miller@gmail.com
# For additional code supporting paper data preparation and pre-analyses see: https://github.com/lark-gorilla/coral_fish 
#
# data (description) included to run code:
# table1_analysis.csv (Summary results from paper Table 1)
# thermal_midpoint_means.csv (average thermal midpoint of tropical species per FG, derived from Stuart Smith et al. 2018)
# JPN_AUS_RMI_CHK_MLD_TMR_trait_master_opt2_clusters.csv (species trait data clustered into functional groups)
# sp_list_summer_only.csv (Summer only Australian species to filter out from trait data)
# Jpn_transects_biomass.csv (transect-level fish biomass data from Japan)
# Aus_transects_biomass.csv (transect-level fish biomass data from Australia)
# Aus_site_species_biomass.csv (site-level fish biomass data from Japan)
# Jpn_site_species_biomass.csv (site-level fish biomass data from Australia)

#rm(list=ls())
#Packages and versions

library(dplyr) #(1.0.6)
library(ggplot2) #(3.3.3 )
library(gridExtra) #(2.3)
library(mice) #(3.13.0)
library(scales) #(1.1.1)
library(ggwordcloud) #(0.5.0) # make sure package dependency "Rccp" is installed and updated
library(vegan) #(2.5-7)
library(ade4) #(1.7-17)
library(emmeans) #(1.6.2-1)
library(ggResidpanel) #(0.3.0)
library(adehabitatHR) #(0.4.19)
library(sf) #(0.9-8)
library(rgeos) #(0.5-5)
library(cluster) #(2.1.2)
library(patchwork) #(1.1.1)

# set working directory
setwd('C:/tropicalisation_paper/NCOMMS-21-26138-T_code_data')

################ read and prepare data ###################
##########################################################

# read in species trait data, with species assigned to 1 of 19 clusters (functional groups;FGs)
dat<-read.csv('JPN_AUS_master_trait.csv', h=T)
dat$Diet<-factor(dat$Diet)
dat$Aggregation<-factor(dat$Aggregation)
dat$Position<-factor(dat$Position)

#Remove Aus summer only species
aus_summer<-read.csv('sp_list_summer_only.csv')
dat[dat$Species %in% aus_summer$Fish,]$AUS_sp<-0

# Read in biomass data
bio_jpn<-read.csv('Jpn_transects_biomass.csv')
bio_aus<-read.csv('Aus_transects_biomass.csv')

# read in species per site biomass data for PCoA
aus_sp_site<-read.csv('Aus_site_species_biomass.csv')
jpn_sp_site<-read.csv('Jpn_site_species_biomass.csv')

# FG to factor
bio_jpn$FG<-factor(bio_jpn$FG)
bio_aus$FG<-factor(bio_aus$FG)

# Remove 2 poorly sampled sites from Jpn and Aus
bio_jpn<-filter(bio_jpn, !SiteID %in% c('JP27', 'JP32'))
bio_aus<-filter(bio_aus, !Site %in% c('Flat Rock', 'Wolf Rock'))

jpn_sp_site<-filter(jpn_sp_site, !SiteID %in% c('JP27', 'JP32'))
aus_sp_site<-filter(aus_sp_site, !Site %in% c('Flat Rock', 'Wolf Rock'))

# correct biomass
bio_aus$cor_biom<-bio_aus$tot_biom/bio_aus$totMsurv
bio_jpn$cor_biom<-bio_jpn$tot_biom/bio_jpn$totMsurv

aus_sp_site$cor_biom<-aus_sp_site$tot_biom/aus_sp_site$totMsurv
jpn_sp_site$cor_biom<-jpn_sp_site$tot_biom/jpn_sp_site$totMsurv

# and summarise by site for sp site data
aus_sp_site<-aus_sp_site%>%group_by(Site, Lat, FG, ThermalAffinity2, Fish)%>%
  summarise(cor_biom=sum(cor_biom))

jpn_sp_site<-jpn_sp_site%>%group_by(SiteID, lat, FG, ThermalAffinity2, SpeciesFish)%>%
  summarise(cor_biom=sum(cor_biom))

########################################


#################### Fig 2 Word Cloud creation ##########################
###################################################################

dat_mice<-mice(dat[c("BodySize","DepthRange", "Diet",  "Position", "Aggregation")], m=5, method=c(rep('norm.predict', 2), rep('polyreg', 3)))
dat_imp<-mice::complete(dat_mice)
dat_imp<-cbind(Species=dat[,1], ThermalAffinity=dat["ThermalAffinity2"], dat_imp, dat[,8:length(dat)])

#manual edit 
dat_imp[dat_imp$Species=='Cantheschenia grandisquamis',]$Diet<-'Omnivore'

# Bodysize

# 6 classes on cuts from Mouillot et al (2014)
dat_imp$BodySize <- cut(dat_imp$BodySize ,breaks = c(-Inf, 7,15,30,50, 80, Inf),
                        labels = c("Tiny","VSmall", "Small", "Medium", "Large", "VLarge"))
# DepthRange

ggplot(data=dat_imp, aes(x=DepthRange))+geom_histogram(binwidth=10, col='black')+
  scale_x_continuous(breaks=seq(0, 420, 10))
# split into those that can only dive to coral reef (~27m) depths, 
# the photic zone (~100m), and deeper > 100m. loosely similar to Mouillot et al (2014) 'Water column' trait
dat_imp$DepthRange <- cut(dat_imp$DepthRange ,breaks = c(-Inf, 30,100, Inf),
                          labels = c("shallow", "mid-depth", "deep"))
dat_imp$FE<-paste(dat_imp$BodySize, dat_imp$DepthRange,
                  dat_imp$Diet, dat_imp$Position, dat_imp$Aggregation)

# How to represent what each FG is using wordclouds of FEs

# reformat data
# split dataframe into list based on rows

FEdat<-dat_imp %>% group_by(FE) %>% summarise(FG=unique(groupk19),num=n())

FEdatlist<-split(FEdat, 1:nrow(FEdat))

FEwordlist<-lapply(FEdatlist, function(x){data.frame(FEcomp=unlist(strsplit(x$FE, ' ')),
                                                     n=x$num, FG=x$FG, FE=x$FE)})
FEword.df<-do.call('rbind', FEwordlist)

FEword.agg<-FEword.df %>% group_by(FG, FEcomp) %>% summarize(sum_word=sum(n)) 

# arranges FEcomp in descending order of n, by FG
FEword.agg<-FEword.agg %>% group_by(FG) %>% arrange(desc(sum_word), .by_group=TRUE)

# help with colouring https://stackoverflow.com/questions/18902485/colored-categories-in-r-wordclouds
FEword.agg$colorlist='blue'
FEword.agg$colorlist<-ifelse(FEword.agg$FEcomp%in%unique(dat_imp$Diet),
                             'orange', FEword.agg$colorlist)
FEword.agg$colorlist<-ifelse(FEword.agg$FEcomp%in%unique(dat_imp$Aggregation),
                             'green', FEword.agg$colorlist)
FEword.agg$colorlist<-ifelse(FEword.agg$FEcomp%in%unique(dat_imp$DepthRange),
                             'purple', FEword.agg$colorlist)
FEword.agg$colorlist<-ifelse(FEword.agg$FEcomp%in%unique(dat_imp$BodySize),
                             'red', FEword.agg$colorlist)

# preserves scaling i.e. sizes of FGs
ggplot(FEword.agg, aes(label=FEcomp, size=sum_word, colour=colorlist))+
  geom_text_wordcloud()+scale_size_area()+facet_wrap(~FG)

FEword.agg.cl<-split(FEword.agg, FEword.agg$FG)

# include seed
out<-lapply(FEword.agg.cl[as.numeric(names(sort(table(dat$groupk19), decreasing = T)))], function(x){
  ggplot(x, aes(label=FEcomp, size=sum_word, colour=colorlist))+
    geom_text_wordcloud(seed=300)+scale_size_area()+theme_minimal()+
    labs(title=paste('FG', unique(x$FG), 'n=', sum(x$sum_word)/5))})

do.call('grid.arrange', out) 

ggsave('C:/coral_fish/outputs/wordclouds.eps',
       plot=do.call('grid.arrange', out),width = 20, height = 20, units = "cm")

# get names for future plots
FEword.agg%>%group_by(FG)%>%
  mutate(totword=sum(sum_word), cumword=cumsum(sum_word), wordperc=cumword/totword*100)%>%
  filter(wordperc<59) ->FGnames

FGnames<-as.data.frame(aggregate(FEcomp~FG, FGnames, 
                                 function(x){paste(x, collapse='-')}))

######################################################################

#### Fig 2 Tropical-subtropical Biomass comparisons (4rt) ####
##############################################################

jpn_cor<-bio_jpn%>%filter(FG %in% c(15, 10, 8, 2, 6, 12, 4, 1, 16))%>%
  dplyr::select(FG, Site.trans.ID,ThermalAffinity2, cor_biom)%>%
  group_by(FG,Site.trans.ID)%>%tidyr::spread(ThermalAffinity2, cor_biom)%>%
  ungroup()%>%group_by(FG)%>%
  summarise(cor_p=cor.test((tropical^0.25), (subtropical^0.25), method = 'kendall')$p.value,
            cop_est=cor.test((tropical^0.25), (subtropical^0.25), method = 'kendall')$estimate)

jpn_cor$txt=ifelse(jpn_cor$cor_p<0.001, paste0('tau=',round(jpn_cor$cop_est, 2), '***'),
                   ifelse(jpn_cor$cor_p<0.01, paste0('tau=',round(jpn_cor$cop_est, 2), '**'),
                          ifelse(jpn_cor$cor_p<0.05, paste0('tau=',round(jpn_cor$cop_est, 2), '*'),
                                 NA)))                     

aus_cor<-bio_aus%>%filter(FG %in% c(15, 10, 8, 2, 6, 12, 4, 1, 16))%>%
  dplyr::select(FG, Site.trans.ID,ThermalAffinity2, cor_biom)%>%
  group_by(FG,Site.trans.ID)%>%tidyr::spread(ThermalAffinity2, cor_biom)%>%
  ungroup()%>%group_by(FG)%>%
  summarise(cor_p=cor.test((tropical^0.25), (subtropical^0.25), method = 'kendall')$p.value,
            cop_est=cor.test((tropical^0.25), (subtropical^0.25), method = 'kendall')$estimate)

aus_cor$txt=ifelse(aus_cor$cor_p<0.001, paste0('tau=',round(aus_cor$cop_est, 2), '***'),
                   ifelse(aus_cor$cor_p<0.01, paste0('tau=',round(aus_cor$cop_est, 2), '**'),
                          ifelse(aus_cor$cor_p<0.05, paste0('tau=',round(aus_cor$cop_est, 2), '*'),
                                 NA)))                     

sg_lat_spans<-data.frame(xmin=c(24.2, 26.2, 28.5, 31,   32.7,  33.38, 34.6), 
                         xmax=c(24.5, 28.4, 30.5, 31.6, 32.82, 33.5, 35 ))


bio_jpn$ThermalAffinity2<-factor(bio_jpn$ThermalAffinity2, levels=c('tropical', 'subtropical'))

jpn_plot_dat<-filter(bio_jpn, FG %in% c(15, 10, 8, 2, 6, 12, 4, 1, 16))

jpn_nsp<-dat%>%filter(JPN_sp==1 & groupk19 %in% c(15, 10, 8, 2, 6, 12, 4, 1, 16))%>%
  group_by(groupk19,ThermalAffinity2)%>%summarise(n())%>%
  dplyr::rename(FG=groupk19, count=`n()`)
jpn_nsp$y=3.5
jpn_nsp$x=34
jpn_nsp[jpn_nsp$ThermalAffinity2=='subtropical',]$x=35

jpn_plot_dat$FG<-factor(jpn_plot_dat$FG, levels=c(15,10,8,2,6,12,4,1,16))
jpn_cor$FG<-factor(jpn_cor$FG, levels=c(15,10,8,2,6,12,4,1,16))
jpn_nsp$FG<-factor(jpn_nsp$FG, levels=c(15,10,8,2,6,12,4,1,16))


p1<-ggplot(data=jpn_plot_dat) + 
  geom_rect(data=sg_lat_spans, aes(ymin=0, ymax=4, xmin=xmin, xmax=xmax), fill='gainsboro')+
  geom_point(data=filter(jpn_plot_dat, ((tot_biom/totMsurv)^0.25)<4),
             aes(x = lat, y = (tot_biom/totMsurv)^0.25, colour=ThermalAffinity2), shape=1)+
  geom_smooth(aes(x = lat, y = (tot_biom/totMsurv)^0.25, colour=ThermalAffinity2),se=F)+
  geom_smooth(data=jpn_plot_dat%>%group_by(Site.trans.ID, FG)%>%summarise(cor_biom=sum(cor_biom), lat=first(lat)),
              aes(x = lat, y = cor_biom^0.25),colour='black', linetype='dotted',se=F)+
  geom_label(data=jpn_cor, aes(x=28.5, y=3.5, label=txt))+
  geom_label(data=jpn_nsp, aes(x=x, y=y, colour=ThermalAffinity2, label=count))+
  facet_grid(FG~., scales='free_y')+
  theme_bw()+theme(legend.position = "none")+
  xlab('Latitude')+ylab('4rt-trans standardised biomass')+
  ggtitle('Japan')+
  theme(strip.background = element_blank(),strip.text.y = element_blank(),
        plot.title = element_text(hjust = 0.5))

sg_lat_spans<-data.frame(xmin=c(-23.4, -24.8, -26.61, -28.19, -29.9, -29.97), 
                         xmax=c(-24.1, -25.3, -26.98, -28.616, -30.96, -30.3))

bio_aus$ThermalAffinity2<-factor(bio_aus$ThermalAffinity2, levels=c('tropical', 'subtropical'))

aus_plot_dat<-filter(bio_aus, FG %in% c(15, 10, 8, 2, 6, 12, 4, 1, 16))

aus_nsp<-dat%>%filter(AUS_sp==1 & groupk19 %in% c(15, 10, 8, 2, 6, 12, 4, 1, 16))%>%
  group_by(groupk19,ThermalAffinity2)%>%summarise(n())%>%
  dplyr::rename(FG=groupk19, count=`n()`)
aus_nsp$y=3.5
aus_nsp$x=-30
aus_nsp[aus_nsp$ThermalAffinity2=='subtropical',]$x=-31

aus_plot_dat$FG<-factor(aus_plot_dat$FG, levels=c(15,10,8,2,6,12,4,1,16))
aus_cor$FG<-factor(aus_cor$FG, levels=c(15,10,8,2,6,12,4,1,16))
aus_nsp$FG<-factor(aus_nsp$FG, levels=c(15,10,8,2,6,12,4,1,16))


p2<-ggplot(data=aus_plot_dat) + 
  geom_rect(data=sg_lat_spans, aes(ymin=0, ymax=4, xmin=xmin, xmax=xmax), fill='gainsboro')+
  geom_point(data=filter(aus_plot_dat,  ((tot_biom/totMsurv)^0.25)<4),
             aes(x = Lat, y = (tot_biom/totMsurv)^0.25, colour=ThermalAffinity2), shape=1)+
  geom_smooth(aes(x = Lat, y = (tot_biom/totMsurv)^0.25, colour=ThermalAffinity2),se=F)+
  geom_smooth(data=aus_plot_dat%>%group_by(Site.trans.ID, FG)%>%summarise(cor_biom=sum(cor_biom), Lat=first(Lat)),
              aes(x = Lat, y = cor_biom^0.25),colour='black', linetype='dotted',se=F)+
  geom_label(data=aus_cor, aes(x=-26, y=3.5, label=txt))+
  geom_label(data=aus_nsp, aes(x=x, y=y, colour=ThermalAffinity2, label=count))+
  facet_grid(FG~., scales='free_y')+scale_x_reverse()+ylab(NULL)+
  ggtitle('Australia')+
  theme_bw()+theme(legend.position = "none")+
  xlab('Latitude')+theme(plot.title = element_text(hjust = 0.5))

#png('C:/coral_fish/outputs/biomass_thermal_plot_sums.png',width = 8, height =12 , units ="in", res =300)
grid.arrange(p2, p1, ncol=2)
dev.off()

#ggsave('C:/coral_fish/outputs/fig2_biomass_labels.eps',
#       plot=grid.arrange(p2, p1, ncol=2),width = 20, height = 29, units = "cm")

################################################################


#### Prepare tropical-only data for tropicalization footprint analyses ####
#########################################################################

# Calc standardisation
# filter out unwanted FGs for per FG objects but not for community total
jpn_trop_prop<-filter(bio_jpn, ThermalAffinity2=='tropical' & lat<25.5)%>%
  group_by(FG)%>%summarise(mean_biom=(mean(cor_biom^0.25, na.rm=T))^4)

aus_trop_prop<-filter(bio_aus, ThermalAffinity2=='tropical' & Lat> -24.5)%>%
  group_by(FG)%>%summarise(mean_biom=(mean(cor_biom^0.25, na.rm=T))^4)

jpn_trop_prop<-left_join(filter(bio_jpn, ThermalAffinity2=='tropical' & 
                                  FG %in% c(15, 10, 8, 2, 6, 12, 4, 1, 16)),
                         jpn_trop_prop, by='FG')

aus_trop_prop<-left_join(filter(bio_aus, ThermalAffinity2=='tropical' & 
                                  FG %in% c(15, 10, 8, 2, 6, 12, 4, 1, 16)),
                         aus_trop_prop, by='FG')

jpn_trop_comm<-filter(bio_jpn, ThermalAffinity2=='tropical')%>%
  group_by(lat, SiteID, Site.trans.ID)%>%summarise(cor_biom=sum(cor_biom)) 
jpn_trop_comm$mean_biom<-as.numeric(filter(jpn_trop_comm, lat<25.5)%>%ungroup()%>%
                                      summarise(mean_biom=(mean(cor_biom^0.25, na.rm=T))^4))

aus_trop_comm<-filter(bio_aus, ThermalAffinity2=='tropical')%>%
  group_by(Lat, Site, Site.trans.ID)%>%summarise(cor_biom=sum(cor_biom))
aus_trop_comm$mean_biom<-as.numeric(filter(aus_trop_comm, Lat> -24.5)%>%ungroup()%>%
                                      summarise(mean_biom=(mean(cor_biom^0.25, na.rm=T))^4))

jpn_trop_prop$FG<-factor(jpn_trop_prop$FG, levels=c(15, 10, 8, 2,6,12,4,1,16))
aus_trop_prop$FG<-factor(aus_trop_prop$FG, levels=c(15, 10, 8, 2,6,12,4,1,16))

# add tropicalization_metric: >1 = more biomass at site compared to tropical site-group
# Using 4th root transformation to make data approx normality (can still handle 0's though)
jpn_trop_comm$trop_met<-(jpn_trop_comm$cor_biom/jpn_trop_comm$mean_biom)^0.25
aus_trop_comm$trop_met<-(aus_trop_comm$cor_biom/aus_trop_comm$mean_biom)^0.25
jpn_trop_prop$trop_met<-(jpn_trop_prop$cor_biom/jpn_trop_prop$mean_biom)^0.25
aus_trop_prop$trop_met<-(aus_trop_prop$cor_biom/aus_trop_prop$mean_biom)^0.25

# create site-groups ("Zones") based on dendrogram clustering 
# Note Zone naming differs slightly from those in paper text
jpn_trop_comm$FG<-'comm'
jpn_trop_prop$FG<-as.character(jpn_trop_prop$FG)
jpn_trop_tests<-rbind(data.frame(jpn_trop_comm), jpn_trop_prop[names(jpn_trop_comm)])
jpn_trop_tests$site.group<-'trop.base'
jpn_trop_tests[jpn_trop_tests$lat > 25 &  jpn_trop_tests$lat < 28.5,]$site.group<-'trop.island'
jpn_trop_tests[jpn_trop_tests$lat > 28.5 &  jpn_trop_tests$lat < 31,]$site.group<-'trans.island'
jpn_trop_tests[jpn_trop_tests$SiteID %in% c('JP28', 'JP29', 'JP30', 'JP31'),]$site.group<-'trans.inland'
jpn_trop_tests[jpn_trop_tests$SiteID %in% c('JP8', 'JP9', 'JP10', 'JP11', 'JP12'),]$site.group<-'trans.headld'
jpn_trop_tests[jpn_trop_tests$lat > 34,]$site.group<-'temp.headld'
table(jpn_trop_tests$SiteID, jpn_trop_tests$site.group)

aus_trop_comm$FG<-'comm'
aus_trop_prop$FG<-as.character(aus_trop_prop$FG)
aus_trop_tests<-rbind(data.frame(aus_trop_comm), aus_trop_prop[names(aus_trop_comm)])
aus_trop_tests$site.group<-'trop.base'
aus_trop_tests[aus_trop_tests$Lat > -25.6 &  aus_trop_tests$Lat < -24.5,]$site.group<-'trans.bay'
aus_trop_tests[aus_trop_tests$Lat > -28 &  aus_trop_tests$Lat < -25.6,]$site.group<-'trans.offshore'
aus_trop_tests[aus_trop_tests$Lat < -28,]$site.group<-'temp.offshore'
aus_trop_tests[aus_trop_tests$Site %in% c('Julian Rock False Trench', 'Julian Rock Nursery', 'Cook Island'),]$site.group<-'trans.temp'
aus_trop_tests[aus_trop_tests$Site %in% c('Muttonbird Island', 'Woolgoolga Reef', 'Woolgoolga Headland', 'North Rock'),]$site.group<-'temp.inshore'

table(aus_trop_tests$Site,aus_trop_tests$site.group)
###################################################################

#### Make FG tropicalization + community trends over latitude (Fig 3 A & C) ####
################################################################################

sg_lat_spans<-data.frame(xmin=c(24.2, 26.2, 28.5, 31,   32.7,  33.38, 34.6), 
                         xmax=c(24.5, 28.4, 30.5, 31.6, 32.82, 33.5, 35 ))

jpn_trop_prop$FG<-factor(jpn_trop_prop$FG, levels=c(15, 10, 8, 2,6,12,4,1,16))
jpn_trop_prop$FG_name<-recode(jpn_trop_prop$FG, 
                              '15'='Benthic Predators', '10'='Upper-benthic Predators',
                              '8'='Benthic Herbivore/Omnivores', '2'='Upper-benthic Planktivores',
                              '6'='Upper-benthic Herbivores','12'='Demersal Predators',
                              '4'='Benthic Planktivores','1'='Upper-benthic Omnivores',
                              '16'='Corallibovores')

jpn_trend<-ggplot(jpn_trop_prop) + 
  geom_rect(data=sg_lat_spans, aes(ymin=0, ymax=2, xmin=xmin, xmax=xmax), fill='grey', alpha=0.5)+
  geom_smooth(aes(x = lat, y = (cor_biom/mean_biom)^0.25, colour=FG_name), se=F)+
  geom_smooth(data=jpn_trop_comm, aes(x = lat, y = (cor_biom/mean_biom)^0.25),
              se=F, colour='black', linetype='dashed')+
  theme_bw()+
  geom_hline(yintercept=0.05^0.25)+scale_x_continuous(breaks=24:35)+
  xlab('Latitude')+ylab('Proportion of tropical biomass')+
  scale_y_continuous(breaks=c(0, 0.05, 0.25, 0.5, 1, 2, 4, 20)^0.25, 
                     labels=c(0,0.05, 0.25, 0.5, 1, 2, 4, 20), minor_breaks = NULL)

sg_lat_spans<-data.frame(xmin=c(-23.4, -24.8, -26.61, -28.19, -29.9, -29.97), 
                         xmax=c(-24.1, -25.3, -26.98, -28.616, -30.96, -30.3))

aus_trop_prop$FG<-factor(aus_trop_prop$FG, levels=c(15, 10, 8, 2,6,12,4,1,16))
aus_trop_prop$FG_name<-recode(aus_trop_prop$FG, 
                              '15'='Benthic Predators', '10'='Upper-benthic Predators',
                              '8'='Benthic Herbivore/Omnivores', '2'='Upper-benthic Planktivores',
                              '6'='Upper-benthic Herbivores','12'='Demersal Predators',
                              '4'='Benthic Planktivores','1'='Upper-benthic Omnivores',
                              '16'='Corallibovores')

aus_trend<-ggplot(aus_trop_prop) + 
  geom_rect(data=sg_lat_spans, aes(ymin=0, ymax=4^0.25, xmin=xmin, xmax=xmax), fill='grey', alpha=0.6)+
  geom_smooth(aes(x = Lat, y = (cor_biom/mean_biom)^0.25, colour=factor(FG)), se=F)+
  geom_smooth(data=aus_trop_comm, aes(x = Lat, y = (cor_biom/mean_biom)^0.25),
              se=F, colour='black', linetype='dashed')+
  theme_bw()+theme(legend.position = "none")+
  geom_hline(yintercept=0.05^0.25)+scale_x_reverse(breaks=-23:-31)+
  xlab('Latitude')+ylab('Proportion of tropical biomass')+
  scale_y_continuous(breaks=c(0, 0.05, 0.25, 0.5, 1, 2, 4)^0.25, 
                     labels=c(0,0.05, 0.25, 0.5, 1, 2, 4), minor_breaks = NULL)
#################################################################################



#### Japan FG tropicalization footprint comparisons (results for Fig 3.B)  ####
###############################################################################

# set comm as intercept
jpn_trop_tests$FG<-factor(jpn_trop_tests$FG, levels=c('comm', 15, 10, 8, 2,6,12,4,1,16))

# remove comm for diff between FG stat reporting
#jpn_trop_tests<-filter(jpn_trop_tests, FG!='comm')
#jpn_trop_tests$FG<-factor(jpn_trop_tests$FG, levels=c( 15, 10, 8, 2,6,12,4,1,16))

## trop.island

ggplot(data=filter(jpn_trop_tests, site.group=='trop.island'), aes(x=FG, y=trop_met))+
  geom_point(aes(colour=SiteID))+geom_boxplot(alpha=0.5) # remember boxplot = medians

m1<-lme(trop_met~FG, random=~1|SiteID, data=filter(jpn_trop_tests, site.group=='trop.island'), 
        weights=varIdent(form=~1|FG))
boxplot(residuals(m1, type='pearson')~filter(jpn_trop_tests, site.group=='trop.island')$FG) # ok good
# pearson and normalized residuals incorporate the effect of the gls variance structure

em1<-emmeans(m1, specs='FG')
pairs(em1)
plot(em1, comparisons = TRUE)

trop_comps_out<-data.frame(site.group='trop.island',data.frame(em1), rbind(c(NA, NA), data.frame(pairs(em1))[1:9, c(1,6)]))

ggplot()+
  geom_jitter(data=filter(jpn_trop_tests, site.group=='trop.island'), aes(x=FG, y=trop_met), shape=1, alpha=0.5, width=0.2, colour='grey')+
  geom_hline(yintercept = trop_comps_out$emmean[1], linetype='dotted')+
  geom_errorbar(data=trop_comps_out, aes(x=FG, ymin=lower.CL, ymax=upper.CL))+
  geom_point(data=trop_comps_out, aes(x=FG, y=emmean, colour=ifelse(p.value>0.05| is.na(p.value), 'blue', 'red')), size=2)+
  theme_bw()+theme(legend.position = "none")

## trans.island

ggplot(data=filter(jpn_trop_tests, site.group=='trans.island'), aes(x=FG, y=trop_met))+
  geom_point(aes(colour=SiteID))+geom_boxplot(alpha=0.5) # remember boxplot = medians

m1<-lme(trop_met~FG, random=~1|SiteID, weights=varIdent(form=~1|FG), data=filter(jpn_trop_tests, site.group=='trans.island'))
#resid_panel(m1)
summary(m1)

em1<-emmeans(m1, specs='FG')
pairs(em1)
plot(em1, comparisons = TRUE)

trop_comps_out<-rbind(trop_comps_out,
                      data.frame(site.group='trans.island',data.frame(em1), rbind(c(NA, NA), data.frame(pairs(em1))[1:9, c(1,6)])))

ggplot()+
  geom_jitter(data=filter(jpn_trop_tests, site.group=='trans.island'), aes(x=FG, y=trop_met), shape=1, alpha=0.5, width=0.2, colour='grey')+
  geom_hline(yintercept =trop_comps_out[trop_comps_out$site.group=='trans.island',]$response[1], linetype='dotted')+
  geom_errorbar(data=filter(trop_comps_out, site.group=='trans.island'), aes(x=FG, ymin=lower.CL, ymax=upper.CL))+
  geom_point(data=filter(trop_comps_out, site.group=='trans.island'), aes(x=FG, y=emmean, colour=ifelse(p.value>0.05| is.na(p.value), 'blue', 'red')), size=2)+
  theme_bw()+theme(legend.position = "none")

## trans.inland

ggplot(data=filter(jpn_trop_tests, site.group=='trans.inland'), aes(x=FG, y=trop_met))+
  geom_point(aes(colour=SiteID))+geom_boxplot(alpha=0.5) # remember boxplot = medians

m1<-lme(trop_met~FG, random=~1|SiteID, weights=varIdent(form=~1|FG), data=filter(jpn_trop_tests, site.group=='trans.inland'))
#resid_panel(m1)
summary(m1)


em1<-emmeans(m1, specs='FG')
pairs(em1)
plot(em1, comparisons = TRUE)

trop_comps_out<-rbind(trop_comps_out,
                      data.frame(site.group='trans.inland',data.frame(em1), rbind(c(NA, NA), data.frame(pairs(em1))[1:9, c(1,6)])))

ggplot()+
  geom_jitter(data=filter(jpn_trop_tests, site.group=='trans.inland'), aes(x=FG, y=trop_met), shape=1, alpha=0.5, width=0.2, colour='grey')+
  geom_hline(yintercept =trop_comps_out[trop_comps_out$site.group=='trans.inland',]$response[1], linetype='dotted')+
  geom_errorbar(data=filter(trop_comps_out, site.group=='trans.inland'), aes(x=FG, ymin=lower.CL, ymax=upper.CL))+
  geom_point(data=filter(trop_comps_out, site.group=='trans.inland'), aes(x=FG, y=emmean, colour=ifelse(p.value>0.05| is.na(p.value), 'blue', 'red')), size=2)+
  theme_bw()+theme(legend.position = "none")


## trans.headld

ggplot(data=filter(jpn_trop_tests, site.group=='trans.headld'), aes(x=FG, y=trop_met))+
  geom_point(aes(colour=SiteID))+geom_boxplot(alpha=0.5) # remember boxplot = medians

m1<-lme(trop_met~FG, random=~1|SiteID, weights=varIdent(form=~1|FG),
        data=filter(jpn_trop_tests, site.group=='trans.headld'),
        control=lmeControl(opt = 'optim')) # different optimizer used to make run
#resid_panel(m1)
summary(m1)


em1<-emmeans(m1, specs='FG')
pairs(em1)
plot(em1, comparisons = TRUE)

trop_comps_out<-rbind(trop_comps_out,
                      data.frame(site.group='trans.headld',data.frame(em1), rbind(c(NA, NA),
                                                                                  data.frame(pairs(em1))[1:9, c(1,6)])))

ggplot()+
  geom_jitter(data=filter(jpn_trop_tests, site.group=='trans.headld'), aes(x=FG, y=trop_met), shape=1, alpha=0.5, width=0.2, colour='grey')+
  geom_hline(yintercept =trop_comps_out[trop_comps_out$site.group=='trans.headld',]$response[1], linetype='dotted')+
  geom_errorbar(data=filter(trop_comps_out, site.group=='trans.headld'), aes(x=FG, ymin=lower.CL, ymax=upper.CL))+
  geom_point(data=filter(trop_comps_out, site.group=='trans.headld'), aes(x=FG, y=emmean, colour=ifelse(p.value>0.05| is.na(p.value), 'blue', 'red')), size=2)+
  theme_bw()+theme(legend.position = "none")

## temp.headld

ggplot(data=filter(jpn_trop_tests, site.group=='temp.headld'), aes(x=FG, y=trop_met))+
  geom_point(aes(colour=SiteID))+geom_boxplot(alpha=0.5) # remember boxplot = medians
# remove FGs with only 0s - these will be sig diff 
zer_fgs<-filter(jpn_trop_tests, site.group=='temp.headld')%>%group_by(FG)%>%summarise(st=sum(trop_met))%>%filter(., st==0)%>%.$FG

m1<-lme(trop_met~FG, random=~1|SiteID, weights=varIdent(form=~1|FG), data=filter(jpn_trop_tests, site.group=='temp.headld' & !FG%in% zer_fgs))
#resid_panel(m1)
summary(m1)

em1<-emmeans(m1, specs='FG')
pairs(em1)
plot(em1, comparisons = TRUE)

trop_comps_out<-rbind(trop_comps_out,
                      data.frame(site.group='temp.headld',data.frame(em1), rbind(c(NA, NA), data.frame(pairs(em1))[1:(9-length(zer_fgs)), c(1,6)])))

ggplot()+
  geom_jitter(data=filter(jpn_trop_tests, site.group=='temp.headld'), aes(x=FG, y=trop_met), shape=1, alpha=0.5, width=0.2, colour='grey')+
  geom_hline(yintercept =trop_comps_out[trop_comps_out$site.group=='temp.headld',]$response[1], linetype='dotted')+
  geom_errorbar(data=filter(trop_comps_out, site.group=='temp.headld'), aes(x=FG, ymin=lower.CL, ymax=upper.CL))+
  geom_point(data=filter(trop_comps_out, site.group=='temp.headld'), aes(x=FG, y=emmean, colour=ifelse(p.value>0.05| is.na(p.value), 'blue', 'red')), size=2)+
  theme_bw()+theme(legend.position = "none")

#### make Fig 3.B to combine with others ####
trop_comps_out<-rbind(trop_comps_out,
                      data.frame(site.group='temp.headld', FG=c(6,12,1,16), emmean=0, SE=0,
                                 df=0, lower.CL=0, upper.CL=0, contrast=NA, p.value=0))

trop_comps_out$site.group<-factor(trop_comps_out$site.group,
                                  levels=c('trop.base', 'trop.island', 'trans.island', 'trans.inland',
                                           'trans.headld', 'temp.headld'))

trop_comps_out$site.group2<-recode(trop_comps_out$site.group,
                                   trop.base = 'Tropical Coral Reef',
                                   trop.island = 'Cold-tropical Reef',
                                   trans.island = 'Transitional Reef',
                                   trans.inland = 'Subtropical Coral Community',
                                   trans.headld = 'Warm-temperate Transitional Reef',
                                   temp.headld = 'Temperate Kelp Reef')



jpn_mods<-ggplot()+
  geom_hline(yintercept=0.05^0.25, linetype='dotted')+
  geom_errorbar(data=trop_comps_out, aes(x=FG, ymin=emmean-SE, ymax=emmean+SE))+
  geom_point(data=trop_comps_out, aes(x=FG, y=emmean, size=ifelse(is.na(p.value), 1, ifelse(p.value>0.05, 1, 2))), shape=1)+
  geom_point(data=trop_comps_out, aes(x=FG, y=emmean, colour=factor(FG)), size=2)+
  theme_bw()+theme(legend.position = "none")+facet_wrap(~site.group2, nrow=1)+
  scale_colour_manual(values = c("black", "#F8766D", "#D39200" ,"#93AA00", "#00BA38",
                                 "#00C19F", "#00B9E3", "#619CFF", "#DB72FB", "#FF61C3"))+
  ylab('Proportion of tropical biomass')+
  scale_y_continuous(breaks=c(0, 0.05, 0.25, 0.5, 1, 2, 4, 20)^0.25, 
                     labels=c(0,0.05, 0.25, 0.5, 1, 2, 4, 20), minor_breaks = NULL)



#### Australia FG tropicalization footprint comparisons (results for Fig 3.D) #### 
##################################################################################

# drop some FGs and set comm as intercept
aus_trop_tests$FG<-factor(aus_trop_tests$FG, levels=c('comm', 15, 10, 8, 2,6,12,4,1,16))

# remove comm for diff between FG stat reporting
#aus_trop_tests<-filter(aus_trop_tests, FG!='comm')
#aus_trop_tests$FG<-factor(aus_trop_tests$FG, levels=c( 15, 10, 8, 2,6,12,4,1,16))


## trans.bay

ggplot(data=filter(aus_trop_tests, site.group=='trans.bay'), aes(x=FG, y=trop_met))+
  geom_point(aes(colour=Site))+geom_boxplot(alpha=0.5) # remember boxplot = medians

# drop 2 big outliers, actually don't
m1<-lme(trop_met~FG, random=~1|Site, weights=varIdent(form=~1|FG), data=filter(aus_trop_tests,
                                                                               site.group=='trans.bay'))
#resid_panel(m1)
summary(m1)

em1<-emmeans(m1, specs='FG')
pairs(em1)

trop_comps_out_aus<-data.frame(site.group='trans.bay',data.frame(em1), rbind(c(NA, NA), data.frame(pairs(em1))[1:9, c(1,6)]))

ggplot()+
  geom_jitter(data=filter(aus_trop_tests, site.group=='trans.bay'), aes(x=FG, y=trop_met), shape=1, alpha=0.5, width=0.2, colour='grey')+
  geom_hline(yintercept =trop_comps_out_aus[trop_comps_out_aus$site.group=='trans.bay',]$response[1], linetype='dotted')+
  geom_errorbar(data=filter(trop_comps_out_aus, site.group=='trans.bay'), aes(x=FG, ymin=lower.CL, ymax=upper.CL))+
  geom_point(data=filter(trop_comps_out_aus, site.group=='trans.bay'), aes(x=FG, y=emmean, colour=ifelse(p.value>0.05| is.na(p.value), 'blue', 'red')), size=2)+
  theme_bw()+theme(legend.position = "none") # stupid massive outlier drags up 12, makes it look ns but it is sig diff


## trans.offshore

ggplot(data=filter(aus_trop_tests, site.group=='trans.offshore'), aes(x=FG, y=trop_met))+
  geom_point(aes(colour=Site))+geom_boxplot(alpha=0.5) # remember boxplot = medians

m1<-lme(trop_met~FG, random=~1|Site, weights=varIdent(form=~1|FG), data=filter(aus_trop_tests, site.group=='trans.offshore'))
#resid_panel(m1)
summary(m1)


em1<-emmeans(m1, specs='FG')
#pairs(em1)
plot(em1, comparisons = TRUE)

trop_comps_out_aus<-rbind(trop_comps_out_aus,
                          data.frame(site.group='trans.offshore',data.frame(em1), rbind(c(NA, NA),
                                                                                        data.frame(pairs(em1))[1:9, c(1,6)])))
ggplot()+
  geom_jitter(data=filter(aus_trop_tests, site.group=='trans.offshore'), aes(x=FG, y=trop_met), shape=1, alpha=0.5, width=0.2, colour='grey')+
  geom_hline(yintercept =trop_comps_out_aus[trop_comps_out_aus$site.group=='trans.offshore',]$response[1], linetype='dotted')+
  geom_errorbar(data=filter(trop_comps_out_aus, site.group=='trans.offshore'), aes(x=FG, ymin=lower.CL, ymax=upper.CL))+
  geom_point(data=filter(trop_comps_out_aus, site.group=='trans.offshore'), aes(x=FG, y=emmean, colour=ifelse(p.value>0.05| is.na(p.value), 'blue', 'red')), size=2)+
  theme_bw()+theme(legend.position = "none")


## trans.temp

ggplot(data=filter(aus_trop_tests, site.group=='trans.temp'), aes(x=FG, y=trop_met))+
  geom_point(aes(colour=Site))+geom_boxplot(alpha=0.5) # remember boxplot = medians

m1<-lme(trop_met~FG, random=~1|Site, weights=varIdent(form=~1|FG), data=filter(aus_trop_tests, site.group=='trans.temp'))
#resid_panel(m1)
summary(m1)

em1<-emmeans(m1, specs='FG')
#pairs(em1)
plot(em1, comparisons = TRUE)

trop_comps_out_aus<-rbind(trop_comps_out_aus,
                          data.frame(site.group='trans.temp',data.frame(em1), rbind(c(NA, NA),
                                                                                    data.frame(pairs(em1))[1:9, c(1,6)])))
ggplot()+
  geom_jitter(data=filter(aus_trop_tests, site.group=='trans.temp'), aes(x=FG, y=trop_met), shape=1, alpha=0.5, width=0.2, colour='grey')+
  geom_hline(yintercept =trop_comps_out_aus[trop_comps_out_aus$site.group=='trans.temp',]$response[1], linetype='dotted')+
  geom_errorbar(data=filter(trop_comps_out_aus, site.group=='trans.temp'), aes(x=FG, ymin=lower.CL, ymax=upper.CL))+
  geom_point(data=filter(trop_comps_out_aus, site.group=='trans.temp'), aes(x=FG, y=emmean, colour=ifelse(p.value>0.05| is.na(p.value), 'blue', 'red')), size=2)+
  theme_bw()+theme(legend.position = "none")

## temp.offshore

ggplot(data=filter(aus_trop_tests, site.group=='temp.offshore'), aes(x=FG, y=trop_met))+
  geom_point(aes(colour=Site))+geom_boxplot(alpha=0.5) # remember boxplot = medians

m1<-lme(trop_met~FG, random=~1|Site, weights=varIdent(form=~1|FG), data=filter(aus_trop_tests, 
                                                                               site.group=='temp.offshore' & trop_met<4)) # remove massive outlier
#resid_panel(m1)
summary(m1)


em1<-emmeans(m1, specs='FG')
#pairs(em1)
plot(em1, comparisons = TRUE)

trop_comps_out_aus<-rbind(trop_comps_out_aus,
                          data.frame(site.group='temp.offshore',data.frame(em1), rbind(c(NA, NA), data.frame(pairs(em1))[1:9, c(1,6)])))

ggplot()+
  geom_jitter(data=filter(aus_trop_tests, site.group=='temp.offshore'), aes(x=FG, y=trop_met), shape=1, alpha=0.5, width=0.2, colour='grey')+
  geom_hline(yintercept =trop_comps_out_aus[trop_comps_out_aus$site.group=='temp.offshore',]$response[1], linetype='dotted')+
  geom_errorbar(data=filter(trop_comps_out_aus, site.group=='temp.offshore'), aes(x=FG, ymin=lower.CL, ymax=upper.CL))+
  geom_point(data=filter(trop_comps_out_aus, site.group=='temp.offshore'), aes(x=FG, y=emmean, colour=ifelse(p.value>0.05| is.na(p.value), 'blue', 'red')), size=2)+
  theme_bw()+theme(legend.position = "none")


## temp.inshore

ggplot(data=filter(aus_trop_tests, site.group=='temp.inshore'), aes(x=FG, y=trop_met))+
  geom_point(aes(colour=Site))+geom_boxplot(alpha=0.5) # remember boxplot = medians
# remove FGs with only 0s - these will be sig diff 
zer_fgs<-filter(aus_trop_tests, site.group=='temp.inshore')%>%group_by(FG)%>%summarise(st=sum(trop_met))%>%filter(., st==0)%>%.$FG

m1<-lme(trop_met~FG, random=~1|Site, weights=varIdent(form=~1|FG), data=filter(aus_trop_tests, site.group=='temp.inshore' & !FG%in% zer_fgs))
#resid_panel(m1)
summary(m1)


em1<-emmeans(m1, specs='FG')
#pairs(em1)
#plot(em1, comparisons = TRUE)

trop_comps_out_aus<-rbind(trop_comps_out_aus,
                          data.frame(site.group='temp.inshore',data.frame(em1), rbind(c(NA, NA), data.frame(pairs(em1))[1:(9-length(zer_fgs)), c(1,6)])))

ggplot()+
  geom_jitter(data=filter(aus_trop_tests, site.group=='temp.inshore'), aes(x=FG, y=trop_met), shape=1, alpha=0.5, width=0.2, colour='grey')+
  geom_hline(yintercept =trop_comps_out_aus[trop_comps_out_aus$site.group=='temp.inshore',]$response[1], linetype='dotted')+
  geom_errorbar(data=filter(trop_comps_out_aus, site.group=='temp.inshore'), aes(x=FG, ymin=lower.CL, ymax=upper.CL))+
  geom_point(data=filter(trop_comps_out_aus, site.group=='temp.inshore'), aes(x=FG, y=emmean, colour=ifelse(p.value>0.05| is.na(p.value), 'blue', 'red')), size=2)+
  theme_bw()+theme(legend.position = "none")

### make Fig 3.D to combine with others ####
trop_comps_out_aus<-rbind(trop_comps_out_aus,
                          data.frame(site.group='temp.inshore', FG=c(6,12,16), emmean=0, SE=0,
                                     df=0, lower.CL=0, upper.CL=0, contrast=NA, p.value=0))


trop_comps_out_aus$site.group<-factor(trop_comps_out_aus$site.group, levels=c('trop.base', 'trans.bay', 'trans.offshore',
                                                                              'trans.temp',  'temp.offshore', 'temp.inshore'))
trop_comps_out_aus$site.group2<-recode(trop_comps_out_aus$site.group,
                                       trop.base = 'Tropical Coral Reef',
                                       trans.bay = 'Cold-tropical Bay',
                                       trans.offshore = 'Warm-transitional Reef',
                                       trans.temp = 'Transitional Reef',
                                       temp.offshore = 'Subtropical Coral Community',
                                       temp.inshore = 'Temperate Kelp Reef')


aus_mods<-ggplot()+
  geom_hline(yintercept=0.05^0.25, linetype='dotted')+
  geom_errorbar(data=trop_comps_out_aus, aes(x=FG, ymin=emmean-SE, ymax=emmean+SE))+
  geom_point(data=trop_comps_out_aus, aes(x=FG, y=emmean, size=ifelse(is.na(p.value), 1, ifelse(p.value>0.05, 1, 2))), shape=1)+
  geom_point(data=trop_comps_out_aus, aes(x=FG, y=emmean, colour=factor(FG)), size=2)+
  theme_bw()+theme(legend.position = "none")+facet_wrap(~site.group2, nrow=1)+
  scale_colour_manual(values = c("black", "#F8766D", "#D39200" ,"#93AA00", "#00BA38",
                                 "#00C19F", "#00B9E3", "#619CFF", "#DB72FB", "#FF61C3"))+
  ylab('Proportion of tropical biomass')+
  scale_y_continuous(breaks=c(0, 0.05, 0.25, 0.5, 1, 2, 4, 20)^0.25, 
                     labels=c(0,0.05, 0.25, 0.5, 1, 2, 4, 20), minor_breaks = NULL)

# write out results
#write.csv(rbind(trop_comps_out, trop_comps_out_aus), 'C:/coral_fish/outputs/sitegroup_FG_tropicalization.csv', quote=F, row.names=F)
############################################################

#### Fig 3 tropicalization footprint plot and write out ####
############################################################


ggsave('C:/coral_fish/outputs/fig3_tropicalization_footprints_legend.eps',
       plot=jpn_trend+jpn_mods+aus_trend+aus_mods+
         plot_layout(ncol=1, nrow=4, guides = 'collect')&theme(legend.position = 'bottom')
       ,width = 21, height = 30, units = "cm")

ggsave('C:/coral_fish/outputs/fig3_tropicalization_footprints_NOlegend.eps',
       plot=jpn_trend+jpn_mods+aus_trend+aus_mods+
         plot_layout(ncol=1, nrow=4, guides = 'collect')&theme(legend.position = 'none')
       ,width = 21, height = 30, units = "cm")

############################################################



#### Look for sig diff between FGs and community-level tropicalization footprint per zone ####
##############################################################################################

jpn_trop_tests$site.group<-factor(jpn_trop_tests$site.group,
                                  levels=c('trop.base', 'trop.island', 'trans.island', 'trans.inland',
                                           'trans.headld', 'temp.headld'))

ggplot(data=filter(jpn_trop_tests, FG=='comm'), aes(x=site.group, y=trop_met))+
  geom_point(aes(colour=SiteID))+geom_boxplot(alpha=0.5) # remember boxplot = medians

m1<-lme(trop_met~site.group, random=~1|SiteID, weights=varIdent(form=~1|site.group),
        data=filter(jpn_trop_tests, FG=='comm'))
resid_panel(m1)
summary(m1)
anova(m1)

em1<-emmeans(m1, specs='site.group')
pairs(em1)
plot(em1, comparisons = TRUE)

aus_trop_tests$site.group<-factor(aus_trop_tests$site.group, levels=c('trop.base', 'trans.bay', 'trans.offshore',
                                                                      'trans.temp',  'temp.offshore', 'temp.inshore'))
ggplot(data=filter(aus_trop_tests, FG=='comm'), aes(x=site.group, y=trop_met))+
  geom_point(aes(colour=Site))+geom_boxplot(alpha=0.5) # remember boxplot = medians

m1<-lme(trop_met~site.group, random=~1|Site, weights=varIdent(form=~1|site.group),
        data=filter(aus_trop_tests, FG=='comm' & !Site.trans.ID%in%
                      c('Julian Rock Nursery_T5','Julian Rock False Trench_T4'))) # removal of 2 big outliers
resid_panel(m1)
summary(m1)
anova(m1)

em1<-emmeans(m1, specs='site.group')
pairs(em1)
plot(em1, comparisons = TRUE)
##############################################################################



#### Test functional groups are independent of thermal preference with #####
################# Stuart-smith (2018) thermal midpoint data ################

# We do not have permission to distribute the Stuart-Smith et al (2018) thermal midpoint data
# however the estimated thermal midpoint per FG, created later in this section 'thermal_midpoint_means.csv',
# is included to allow interrogation of results.

therm_mid<-readxl::read_xlsx('C:/coral_fish/sourced_data/stuart_smith_thermal_midpoints/Thermal niche midpoints.xlsx')
tanika_match<-read.csv('C:/coral_fish/sourced_data/stuart_smith_thermal_midpoints/species_match_tanika.csv')
# bodge loop to fix 13 species names
for(i in tanika_match[tanika_match$Name.in.st.data!='',]$Name.in.st.data){
  therm_mid[therm_mid$SPECIES_NAME==i,]$SPECIES_NAME<-as.character(tanika_match[tanika_match$Name.in.st.data==i,]$Not.matched)}

dat_ss<-left_join(dat, therm_mid, by=c('Species'='SPECIES_NAME'))
which(is.na(dat_ss$`95th SSTmax`)) # some naming mis-matches/ missing sp ~80
table(dat_ss$ThermalAffinity2, dat_ss$`Temp-Trop (23cutoff)`) # some differences

# check for coverage of FGs and compare within each region
# filter to FGs we're interested in
dat_ss<-filter(dat_ss, groupk19 %in% c(15, 10, 8, 2, 6, 12, 4, 1, 16))
dat_ss$groupk19<-factor(dat_ss$groupk19)
dat_ss$sst95<-dat_ss$`95th SSTmax`

#JPN
dat_ss%>%filter(JPN_sp>0 & ThermalAffinity2=='tropical')%>%group_by(groupk19)%>%
  summarise(n_sp=n(), n_sp_conf=length(confidence[which(confidence>1)])) # pretty good

conf_tm_jpn<-dat_ss%>%filter(JPN_sp>0 & ThermalAffinity2=='tropical' & confidence>1)

qplot(data=conf_tm_jpn, x=groupk19, y=sst95, geom='boxplot') # looks good but need closer 

m1<-lm(sst95~groupk19, data=conf_tm_jpn[-80,]) #remove outlier
#par(mfrow=c(2,2));plot(m1)
#resid_panel(m1)
#m1.1<-lm(sst95~groupk19, data=conf_tm_jpn[-c(80, which(resid(m1)>1|resid(m1)< -1)),])
resid_panel(m1.1)
summary(m1)
anova(m1) # ns
pairs(emmeans(m1, 'groupk19'))
jpn_mean<-emmeans(m1, 'groupk19')
boxplot(resid(m1, type='pearson')~factor(conf_tm_jpn$groupk19))

ggplot(data=data.frame(emmeans(m1, 'groupk19')), aes(x=factor(groupk19), y=emmean))+geom_pointrange(aes(ymin=lower.CL, ymax=upper.CL))+xlab('Functional Group')+ylab('Thermal Midpoint')

#try with GLS to be sure
m2<-gls(sst95~factor(groupk19), data=conf_tm_jpn,
        weights=varIdent(form=~1|groupk19))
boxplot(resid(m2, type='pearson')~factor(conf_tm_jpn$groupk19))
# no better, remember these groups have different n() so standard 
# error different anyway
p1<-ggplot(data=data.frame(emmeans(m1, 'groupk19')), aes(x=factor(groupk19), y=emmean))+
  geom_pointrange(aes(ymin=lower.CL, ymax=upper.CL))+xlab('Functional Group')+
  ylab('Thermal Midpoint (°C)')+scale_y_continuous(limits=c(30.5, 31.6),
                                                   breaks=c(30.5,30.75, 31, 31.25, 31.5, 31.5))+labs(title ='Japan')

#AUS

dat_ss%>%filter(AUS_sp>0 & ThermalAffinity2=='tropical')%>%group_by(groupk19)%>%
  summarise(n_sp=n(), n_sp_conf=length(confidence[which(confidence>1)])) # pretty good

conf_tm_aus<-dat_ss%>%filter(AUS_sp>0 & ThermalAffinity2=='tropical' & confidence>1)

qplot(data=conf_tm_aus, x=groupk19, y=sst95, geom='boxplot') # looks good but need closer 

m1<-lm(sst95~groupk19, data=conf_tm_aus[-c(18, 186, 89),]) #remove outliers
par(mfrow=c(2,2));plot(m1)
resid_panel(m1)
summary(m1)
anova(m1) # ns
pairs(emmeans(m1, 'groupk19'))
aus_mean<-emmeans(m1, 'groupk19')
boxplot(resid(m1, type='pearson')~factor(conf_tm_aus$groupk19))

p2<-ggplot(data=data.frame(emmeans(m1, 'groupk19')), aes(x=factor(groupk19), y=emmean))+
  geom_pointrange(aes(ymin=lower.CL, ymax=upper.CL))+xlab('Functional Group')+
  ylab('Thermal Midpoint (°C)')+scale_y_continuous(limits=c(30.5, 31.6),
                                                   breaks=c(30.5,30.75, 31, 31.25, 31.5, 31.5))+labs(title ='Australia')

#png('C:/coral_fish/outputs/themal_midpoint_suppl.png',width =8, height =4 , units ="in", res =600)
grid.arrange(p1, p2, ncol=2)
dev.off()

write.csv(rbind(data.frame(region='Japan', jpn_mean),data.frame(region='Australia', aus_mean)),
          'C:/coral_fish/outputs/thermal_midpoint_means.csv', quote=F, row.names=F)

## Check for site.group differences between FG that could be explained by sst differences

jpn_sp_site$site.group<-'trop.base'
jpn_sp_site[jpn_sp_site$lat > 25 &  jpn_sp_site$lat < 28.5,]$site.group<-'trop.island'
jpn_sp_site[jpn_sp_site$lat > 28.5 &  jpn_sp_site$lat < 31,]$site.group<-'trans.island'
jpn_sp_site[jpn_sp_site$SiteID %in% c('JP28', 'JP29', 'JP30', 'JP31'),]$site.group<-'trans.inland'
jpn_sp_site[jpn_sp_site$SiteID %in% c('JP8', 'JP9', 'JP10', 'JP11', 'JP12'),]$site.group<-'trans.headld'
jpn_sp_site[jpn_sp_site$lat > 34,]$site.group<-'temp.headld'

aus_sp_site$site.group<-'trop.base'
aus_sp_site[aus_sp_site$Lat > -25.6 &  aus_sp_site$Lat < -24.5,]$site.group<-'trans.bay'
aus_sp_site[aus_sp_site$Lat > -28 &  aus_sp_site$Lat < -25.6,]$site.group<-'trans.offshore'
aus_sp_site[aus_sp_site$Lat < -28,]$site.group<-'temp.offshore'
aus_sp_site[aus_sp_site$Site %in% c('Julian Rock False Trench', 'Julian Rock Nursery', 'Cook Island'),]$site.group<-'trans.temp'
aus_sp_site[aus_sp_site$Site %in% c('Muttonbird Island', 'Woolgoolga Reef', 'Woolgoolga Headland', 'North Rock'),]$site.group<-'temp.inshore'

aus_sp_site$site.group<-factor(aus_sp_site$site.group,
                               levels=c('trop.base', 'trans.bay', 'trans.offshore','trans.temp',  'temp.offshore', 'temp.inshore'))
jpn_sp_site$site.group<-factor(jpn_sp_site$site.group,
                               levels=c('trop.base', 'trop.island', 'trans.island', 'trans.inland','trans.headld', 'temp.headld'))

jpn_sp_site<-jpn_sp_site%>%filter(FG %in% c(15, 10, 8, 2, 6, 12, 4, 1, 16)
                                  & ThermalAffinity2=='tropical')
aus_sp_site<-aus_sp_site%>%filter(FG %in% c(15, 10, 8, 2, 6, 12, 4, 1, 16)
                                  & ThermalAffinity2=='tropical')

jpn_sp_site<-left_join(jpn_sp_site, dat_ss%>%dplyr::select(Species, confidence, sst95), 
                       by=c('SpeciesFish'='Species'))
aus_sp_site<-left_join(aus_sp_site, dat_ss%>%dplyr::select(Species, confidence, sst95), 
                       by=c('Fish'='Species'))

ggplot(data=filter(jpn_sp_site, confidence>1))+
  geom_boxplot(aes(y=sst95, x=factor(FG)))+facet_wrap(~site.group, scales='free')

ggplot(data=filter(jpn_sp_site, confidence>1))+
  geom_boxplot(aes(y=sst95, x=factor(site.group)))+facet_wrap(~FG, scales='free_y')

em1<-filter(jpn_sp_site, confidence>1) %>% group_by(site.group) %>% do(tidy(emmeans(lm(sst95 ~ factor(FG), .), 'FG')))

filter(jpn_sp_site, confidence>1) %>% group_by(site.group) %>% do(tidy(anova(lm(sst95 ~ factor(FG), .) )))
#trop.base = 0.018, trans.island=0.013

p1<-ggplot()+geom_jitter(data=filter(jpn_sp_site, confidence>1), 
                         aes(x=factor(FG),y=sst95), shape=1, alpha=0.5, colour='red', width=0.1)+
  geom_pointrange(data=em1, aes(x=factor(FG),y=estimate, ymin=conf.low, ymax=conf.high))+
  facet_wrap(~site.group, scales='free', nrow=1)+theme_bw()+
  xlab('Functional Group')+ylab('Thermal Midpoint (°C)')+ labs(title ='Japan')


ggplot(data=filter(aus_sp_site, confidence>1))+
  geom_boxplot(aes(y=sst95, x=factor(FG)))+facet_wrap(~site.group, scales='free')

ggplot(data=filter(aus_sp_site, confidence>1))+
  geom_boxplot(aes(y=sst95, x=factor(site.group)))+facet_wrap(~FG, scales='free_y')

em1<-filter(aus_sp_site, confidence>1) %>% group_by(site.group) %>% do(tidy(emmeans(lm(sst95 ~ factor(FG), .), 'FG')))

filter(aus_sp_site, confidence>1) %>% group_by(site.group) %>% do(tidy(anova(lm(sst95 ~ factor(FG), .) )))
# trop.base=0.005, trans.bay <0.001, temp.inshore <0.001

p2<-ggplot()+geom_jitter(data=filter(aus_sp_site, confidence>1), 
                         aes(x=factor(FG),y=sst95), shape=1, alpha=0.5, colour='red', width=0.1)+
  geom_pointrange(data=em1, aes(x=factor(FG),y=estimate, ymin=conf.low, ymax=conf.high))+
  facet_wrap(~site.group, scales='free', nrow=1)+theme_bw()+
  xlab('Functional Group')+ylab('Thermal Midpoint (°C)')+ labs(title ='Australia')

#png('C:/coral_fish/outputs/themal_midpoint_zone_suppl.png',width =8, height =8 , units ="in", res =600)
grid.arrange(p1, p2, nrow=2)
dev.off()
###########################################################################

#### Trait space setup and data preparation ######
##################################################

# recreate distance matrix from clustering

#order aggregation
dat$Aggregation<-factor(dat$Aggregation, levels=c("solitary", "pairs","groups","schools"), ordered = T)

distlog<-daisy(dat[,c("BodySize","Diet",  "Position", "Aggregation", 'DepthRange')],
               metric = "gower",stand = FALSE, type = list(logratio = c(1,5)))

func_dudi<-dudi.pco(d = cailliez(distlog, print=TRUE, cor.zero = F), scannf = FALSE, nf = 4)
#func_dudi<-dudi.pco(d = (distlog+3.04421), scannf = FALSE, nf =4)

screeplot(func_dudi)
hist(func_dudi$eig)

func_pco<-data.frame(func_dudi$li,dat)

#### Quality of dendrogram ######
#################################
# Variance explained contribution of PCoA first 40 axes
100* func_dudi$eig[1:40]/sum(func_dudi$eig[func_dudi$eig>0.007])


library(mFD)

d1<-dat
d1$BodySize<-log(dat$BodySize)
d1$DepthRange<-log(dat$DepthRange)
distlog1<-daisy(d1[,c("BodySize","Diet",  "Position", "Aggregation", 'DepthRange')], metric = "gower",stand = FALSE)

d1_cat<-data.frame(trait_name =c("BodySize","Diet",  "Position", "Aggregation", 'DepthRange'),
                   trait_type=c('Q','N', 'N', 'O', 'Q' ))
row.names(d1)<-d1$Species
sp_dist_fish <- mFD::funct.dist(
  sp_tr         = d1[,c("BodySize","Diet",  "Position", "Aggregation", 'DepthRange')],
  tr_cat        = d1_cat,
  metric        = "gower", stop_if_NA=F)

fish_quality <- mFD::quality.fspaces(
  sp_dist             = sp_dist_fish,
  maxdim_pcoa         = 2,
  deviation_weighting = "absolute",
  fdist_scaling       = F,
  fdendro             = "average")

round(fish_quality$"quality_fspaces", 3)   

# plot it
btree<-hclust(distlog, method='average')
copo<-cophenetic(btree)

dl1 <- lower.tri(distlog, diag=TRUE)
cp1 <- lower.tri(copo, diag=TRUE)
df1<-data.frame(gower_d=distlog[dl1],
           cophenetic_d=copo[cp1])

ggplot(data=df1)+geom_point(aes(x=gower_d, y=cophenetic_d))

df2<-data.frame(m1=rep(1:nrow(dat), nrow(dat)), m2=sort(rep(1:nrow(dat), nrow(dat))),
           gower_d=as.vector(as.matrix(distlog)), copo_d=as.vector(as.matrix(copo))) 

dtemp<-data.frame(id=as.integer(row.names(dat)), groupk19=dat$groupk19, Species=dat$Species)
df3<-left_join(df2, dtemp, by=c("m1"="id"))
df4<-left_join(df3, dtemp, by=c("m2"="id"))

df5<-df4%>%filter(groupk19.x %in% c(15, 10, 8, 2, 6, 12, 4, 1, 16))%>%
  filter(groupk19.y %in% c(15, 10, 8, 2, 6, 12, 4, 1, 16))

#write.csv(df5, "C:/tropicalisation_paper/submission_4th_round/gower_cophenetic_pairs.csv", quote=F, row.names=F)

# show all sp pairs with Gower < 0.2 and that are in different groups - all due to one species "Cantheschenia grandisquamis"
df6<-df5%>%filter(gower_d<0.2)%>%filter(groupk19.x!=groupk19.y)
df6a<-df5%>%filter(gower_d<0.2)%>%filter(groupk19.x==groupk19.y)

df7<-df6%>%filter(Species.x=="Cantheschenia grandisquamis")
df7$Pair=paste(df7$Species.x, df7$Species.y, sep="-")
df7$Group.Pair=paste(df7$groupk19.x, df7$groupk19.y, sep="-")

library(ggpointdensity)
library(viridis)
ggplot()+
  geom_abline(intercept = 0, slope = 1)+
  geom_pointdensity(data=df5, mapping = aes(x = gower_d, y = copo_d),show.legend = F) +
  geom_point(data=df6a, aes(x = gower_d, y = copo_d),colour="black", shape=1)+
  geom_point(data=df7, aes(x = gower_d, y = copo_d, colour=Group.Pair))
  scale_color_viridis()+xlab("Gower distance")+ylab("Cophenetic distance")+theme_bw()

#################################

# see how it looks

ppp<- ggplot()+
  geom_hline(yintercept=0, linetype="dotted") + 
  geom_vline(xintercept=0,  linetype="dotted")

ppp+geom_point(data=func_pco, aes(x=A1, y=A2))

# make plots per site-groups

jpn_sp_site$site.group<-'trop.base'
jpn_sp_site[jpn_sp_site$lat > 25 &  jpn_sp_site$lat < 28.5,]$site.group<-'trop.island'
jpn_sp_site[jpn_sp_site$lat > 28.5 &  jpn_sp_site$lat < 31,]$site.group<-'trans.island'
jpn_sp_site[jpn_sp_site$SiteID %in% c('JP28', 'JP29', 'JP30', 'JP31'),]$site.group<-'trans.inland'
jpn_sp_site[jpn_sp_site$SiteID %in% c('JP8', 'JP9', 'JP10', 'JP11', 'JP12'),]$site.group<-'trans.headld'
jpn_sp_site[jpn_sp_site$lat > 34,]$site.group<-'temp.headld'
table(jpn_sp_site$SiteID, jpn_sp_site$site.group)

aus_sp_site$site.group<-'trop.base'
aus_sp_site[aus_sp_site$Lat > -25.6 &  aus_sp_site$Lat < -24.5,]$site.group<-'trans.bay'
aus_sp_site[aus_sp_site$Lat > -28 &  aus_sp_site$Lat < -25.6,]$site.group<-'trans.offshore'
aus_sp_site[aus_sp_site$Lat < -28,]$site.group<-'temp.offshore'
aus_sp_site[aus_sp_site$Site %in% c('Julian Rock False Trench', 'Julian Rock Nursery', 'Cook Island'),]$site.group<-'trans.temp'
aus_sp_site[aus_sp_site$Site %in% c('Muttonbird Island', 'Woolgoolga Reef', 'Woolgoolga Headland', 'North Rock'),]$site.group<-'temp.inshore'

# join with PCoA axis data

# leave biomass summed per site

jpn_sp_site_pco<-left_join(jpn_sp_site, func_pco[,1:5], by=c('SpeciesFish'='Species'))
aus_sp_site_pco<-left_join(aus_sp_site, func_pco[,1:5], by=c('Fish'='Species'))

names(jpn_sp_site_pco)[1]<-'Site'
names(jpn_sp_site_pco)[2]<-'Lat'
######################################################################################

##### Custom function to make kernels in trait space surrounding species points #######
##### grouped by FG and thermal affinity, then calculate overlap within each Zone #####
#######################################################################################
funcOvl<-function(pcodat=mydat, FGz=c(10, 15), bywhat='site', mkern=F, selPC1="A5", selPC2="A6")
{  pcodat<-as.data.frame(pcodat)

if(bywhat=='site'){
  pcodat$ID<-paste(pcodat$Site, pcodat$FG, pcodat$ThermalAffinity2, sep='@')}
if(bywhat=='site.group'){
  pcodat$ID<-paste(pcodat$site.group, pcodat$FG, pcodat$ThermalAffinity2, sep='@')}

# loop to find IDs with < 5 data points and replicate until n=5 so kernels can be made
for(i in data.frame(table(pcodat$ID))[which(data.frame(table(pcodat$ID))$Freq<5),]$Var1)
{
  pcodat<-rbind(pcodat, 
                do.call(rbind,replicate(4, pcodat[pcodat$ID==i,], simplify=F)))
}

site_ovl_comp<-NULL
site_ovl_filt<-NULL
site_ovl_kernz<-NULL
for(j in FGz)
{
  pcFG<-pcodat[pcodat$FG==j,]
  
  if(bywhat=='site'){
    temp_base<-pcFG[pcFG$site.group=='trop.base',]
    temp_base$ID<-paste(temp_base$site.group, temp_base$FG, temp_base$ThermalAffinity2, sep='@')
    pcFG<-rbind(pcFG, temp_base)}  
  
  spdf<-SpatialPointsDataFrame(coords=pcFG[c(selPC1,selPC2)],  data=data.frame(ID=pcFG$ID))
  
  spdf$ID<-factor(spdf$ID)
  
  KDE.Surface <- kernelUD(spdf,same4all = T, h=0.075, grid=500)
  
  if(mkern==T){
    KDE.99_site <- getverticeshr(KDE.Surface, percent = 99)
    KDE.99_site<-st_as_sf(KDE.99_site)
    temp<- do.call(rbind, lapply(strsplit(as.character(KDE.99_site$id), '@'),
                                 function(x) data.frame(Site=x[1], FG=x[2], ThermalAffinity2=x[3])))
    KDE.99_site$Site<-temp$Site
    KDE.99_site$FG<-temp$FG
    KDE.99_site$ThermalAffinity2<-temp$ThermalAffinity2
  }else{KDE.99_site<-NULL}
  
  ovl<-kerneloverlaphr(KDE.Surface, percent=99, meth="HR")
  ovl<-as.data.frame.table(ovl)
  ovl<-ovl[-which(ovl$Var1==ovl$Var2),]
  ovl<-cbind(ovl,do.call(rbind, lapply(strsplit(as.character(ovl$Var1), '@'),
                                       function(x) data.frame(V1Site=x[1], V1FG=x[2], V1ThermalAffinity2=x[3]))))
  ovl<-cbind(ovl,do.call(rbind, lapply(strsplit(as.character(ovl$Var2), '@'),
                                       function(x) data.frame(V2Site=x[1], V2FG=x[2], V2ThermalAffinity2=x[3]))))
  ovl<-ovl[ovl$V1FG==ovl$V2FG,]
  
  ovl_comp<-ovl[ovl$V1Site==ovl$V2Site,]
  ovl_comp<-ovl_comp%>%group_by(V1Site, V1FG)%>%summarise_all(first)%>%as.data.frame() # prop of tropical FG space covered by subtropical FG
  names(ovl_comp)[names(ovl_comp)=="V1Site"]<-'Site'
  names(ovl_comp)[names(ovl_comp)=="V1FG"]<-'FG'
  if(bywhat=='site'){ovl_comp<-ovl_comp[ovl_comp$Site!='trop.base',]}
  
  ovl_filt<-ovl[ovl$V1ThermalAffinity2==ovl$V2ThermalAffinity2,]
  ovl_filt<-ovl_filt[ovl_filt$V1ThermalAffinity2=='tropical',]
  ovl_filt<-ovl_filt[ovl_filt$V1Site=='trop.base',]
  names(ovl_filt)[names(ovl_filt)=="V2Site"]<-'Site'
  names(ovl_filt)[names(ovl_filt)=="V1FG"]<-'FG'
  
  site_ovl_comp<-rbind(site_ovl_comp, ovl_comp)
  site_ovl_filt<-rbind(site_ovl_filt, ovl_filt)
  site_ovl_kernz<-rbind(site_ovl_kernz, KDE.99_site)
  print(j)
}
return(list(site_ovl_comp, site_ovl_filt, site_ovl_kernz))
} # function end #

# run function per site for lat plots with tropicalization metric
ovl_site_aus<-funcOvl(pcodat=aus_sp_site_pco, FGz=c(15, 10, 8, 2,6, 12, 4,1, 16), bywhat='site', mkern=F, selPC1="A1", selPC2="A2")
ovl_site_jpn<-funcOvl(pcodat=jpn_sp_site_pco, FGz=c(15, 10, 8, 2,6, 12, 4,1, 16), bywhat='site', mkern=F, selPC1="A1", selPC2="A2")

# run function per site.group for kernel visualisation in FG space
ovl_site.group_aus<-funcOvl(pcodat=aus_sp_site_pco, FGz=c(15, 10, 8, 2,6, 12, 4,1, 16), bywhat='site.group', mkern=T, selPC1="A1", selPC2="A2")
ovl_site.group_jpn<-funcOvl(pcodat=jpn_sp_site_pco, FGz=c(15, 10, 8, 2,6, 12, 4,1, 16), bywhat='site.group', mkern=T, selPC1="A1", selPC2="A2")

# fill overlap with 0's where no overlap
ovl_site_aus[[1]]<-ovl_site_aus[[1]]%>%tidyr::complete(Site, FG, fill=list(Freq=0))
ovl_site_aus[[2]]<-ovl_site_aus[[2]]%>%tidyr::complete(Site, FG, fill=list(Freq=0))

ovl_site_jpn[[1]]<-ovl_site_jpn[[1]]%>%tidyr::complete(Site, FG, fill=list(Freq=0))
ovl_site_jpn[[2]]<-ovl_site_jpn[[2]]%>%tidyr::complete(Site, FG, fill=list(Freq=0))
########################################################################################


#### Create plots (Fig 4) of functional niche overlap and change (competition and filtering) ####
#################################################################################################

## AUSTRALIA
aus_ovl_comp<-left_join(ovl_site_aus[[1]],aus_sp_site%>%group_by(Site)%>%
                          summarise_all(first)%>%dplyr::select(Site, Lat), by='Site')
aus_ovl_filt<-left_join(ovl_site_aus[[2]],aus_sp_site%>%group_by(Site)%>%
                          summarise_all(first)%>%dplyr::select(Site, Lat), by='Site')

aus_ovl_comp<-aus_ovl_comp[-which(is.na(aus_ovl_comp$Site)),] # remove na row from FG 16 no overlap

sg_lat_spans<-data.frame(xmin=c(-23.4, -24.8, -26.61, -28.19, -29.9, -29.97), 
                         xmax=c(-24.1, -25.3, -26.98, -28.616, -30.96, -30.3))
p1<-ggplot()+
  geom_rect(data=sg_lat_spans, aes(ymin=0, ymax=2, xmin=xmin, xmax=xmax), fill='grey', alpha=0.5)+
  geom_smooth(data=aus_trop_prop, aes(x=Lat, y=trop_met), se=F, colour='black')+
  geom_smooth(data=aus_ovl_comp, aes(x=Lat, y=Freq), colour='red', se=F)+
  geom_smooth(data=aus_ovl_filt, aes(x=Lat, y=Freq), colour='blue', se=F)+
  geom_point(data=filter(aus_trop_prop, trop_met<2), aes(x=Lat, y=trop_met))+
  geom_point(data=aus_ovl_comp, aes(x=Lat, y=Freq), colour='red')+
  geom_point(data=aus_ovl_filt, aes(x=Lat, y=Freq), colour='blue')+
  theme_bw()+facet_grid(FG~.)+scale_x_reverse()
# note removal of outliers on geom_point: allows smooth tofit to full data
# but not show these outliers on plot

levz<-c('trop.base', 'trans.bay', 'trans.offshore', 'trans.temp',
        'temp.inshore', 'temp.offshore')

aus_sp_site_pco$site.group<-factor(aus_sp_site_pco$site.group,levels=levz)
names(ovl_site.group_aus[[3]])[names(ovl_site.group_aus[[3]])=='Site']<-'site.group'
names(ovl_site.group_aus[[2]])[names(ovl_site.group_aus[[2]])=='Site']<-'site.group'
names(ovl_site.group_aus[[1]])[names(ovl_site.group_aus[[1]])=='Site']<-'site.group'

ovl_site.group_aus[[3]]$site.group<-factor(ovl_site.group_aus[[3]]$site.group,levels=levz)
ovl_site.group_aus[[1]]$site.group<-factor(ovl_site.group_aus[[1]]$site.group,levels=levz)
ovl_site.group_aus[[2]]$site.group<-factor(ovl_site.group_aus[[2]]$site.group,levels=levz)

#ovl_site.group_aus[[1]]<-ovl_site.group_aus[[1]][-which(is.na(ovl_site.group_aus[[1]]$site.group)),]

ap1<-aus_sp_site_pco[aus_sp_site_pco$FG %in% c(15, 10, 8, 2,6, 12, 4,1, 16),]
ap1$FG<-factor(ap1$FG)

ap1$ThermalAffinity2<-factor(ap1$ThermalAffinity2, levels=c('tropical', 'subtropical'))
ovl_site.group_aus[[3]]$ThermalAffinity2<-factor(ovl_site.group_aus[[3]]$ThermalAffinity2, levels=c('tropical', 'subtropical'))

p2<-ppp+geom_point(data=ap1, aes(x=A1, y=A2, colour=ThermalAffinity2))+
  geom_sf(data=ovl_site.group_aus[[3]], aes(fill=ThermalAffinity2), alpha=0.5)+
  geom_text(data=ovl_site.group_aus[[1]], aes(x=0.9, y=0.9, label=paste(round(Freq, 2)*100,'%', sep='')))+
  geom_text(data=ovl_site.group_aus[[2]], aes(x=0.7, y=0.7, label=paste(round(Freq, 2)*100,'%', sep='')), colour='red')+
  scale_y_continuous(paste('PC2', sprintf('(%0.1f%% explained var.)',
                                          100* func_dudi$eig[2]/sum(func_dudi$eig[func_dudi$eig>0.007]))))+
  scale_x_continuous(paste('PC1', sprintf('(%0.1f%% explained var.)',
                                          100* func_dudi$eig[1]/sum(func_dudi$eig[func_dudi$eig>0.007]))))+
  facet_grid(FG~site.group)+theme_minimal()+theme(legend.position = "none")

p2a<-ppp+
  geom_sf(data=ovl_site.group_aus[[3]], aes(colour=ThermalAffinity2), fill=NA)+
  geom_text(data=ovl_site.group_aus[[1]], aes(x=0.9, y=0.9, label=paste(round(Freq, 2)*100,'%', sep='')))+
  geom_text(data=ovl_site.group_aus[[2]], aes(x=0.7, y=0.7, label=paste(round(Freq, 2)*100,'%', sep='')), colour='red')+
  scale_y_continuous(paste('PC2', sprintf('(%0.1f%% explained var.)',
                                          100* func_dudi$eig[2]/sum(func_dudi$eig[func_dudi$eig>0.007]))))+
  scale_x_continuous(paste('PC1', sprintf('(%0.1f%% explained var.)',
                                          100* func_dudi$eig[1]/sum(func_dudi$eig[func_dudi$eig>0.007]))))+
  facet_grid(FG~site.group)+theme_minimal()+theme(legend.position = "none")


## JAPAN

jpn_ovl_comp<-left_join(ovl_site_jpn[[1]],jpn_sp_site%>%group_by(SiteID)%>%summarise_all(first)%>%dplyr::select(SiteID, lat), by=c('Site'='SiteID'))
jpn_ovl_filt<-left_join(ovl_site_jpn[[2]],jpn_sp_site%>%group_by(SiteID)%>%summarise_all(first)%>%dplyr::select(SiteID, lat), by=c('Site'='SiteID'))

# THESE NOT IN AUS
jpn_trop_prop$FG<-factor(jpn_trop_prop$FG,
                         levels=c(15, 10, 8, 2,6, 12, 4,1, 16))
jpn_ovl_comp$FG<-factor(jpn_ovl_comp$FG,
                        levels=c(15, 10, 8, 2,6, 12, 4,1, 16))
jpn_ovl_filt$FG<-factor(jpn_ovl_filt$FG,
                        levels=c(15, 10, 8, 2,6, 12, 4,1, 16))

jpn_ovl_comp<-jpn_ovl_comp[-which(is.na(jpn_ovl_comp$Site)),]

sg_lat_spans<-data.frame(xmin=c(24.2, 26.2, 28.5, 31,   32.7,  33.38, 34.6), 
                         xmax=c(24.5, 28.4, 30.5, 31.6, 32.82, 33.5, 35 ))

p3<-ggplot()+
  geom_rect(data=sg_lat_spans, aes(ymin=0, ymax=2, xmin=xmin, xmax=xmax), fill='grey', alpha=0.5)+
  geom_smooth(data=jpn_trop_prop, aes(x=lat, y=trop_met), se=F, colour='black')+
  geom_smooth(data=jpn_ovl_comp, aes(x=lat, y=Freq), colour='red', se=F)+
  geom_smooth(data=jpn_ovl_filt, aes(x=lat, y=Freq), colour='blue', se=F)+
  geom_point(data=jpn_trop_prop[jpn_trop_prop$trop_met<2,], aes(x=lat, y=trop_met))+
  geom_point(data=jpn_ovl_comp, aes(x=lat, y=Freq), colour='red')+
  geom_point(data=jpn_ovl_filt, aes(x=lat, y=Freq), colour='blue')+
  theme_bw()+facet_grid(FG~.)
# FYI tropicalization curve is same when using site-aggregated mean data
# jpn_trop_prop%>%group_by(FG, SiteID)%>%summarise_all(mean)

levz<-c('trop.base', 'trop.island', 'trans.island', 'trans.inland',
        'trans.headld', 'temp.headld')

jpn_sp_site_pco$site.group<-factor(jpn_sp_site_pco$site.group,levels=levz)
names(ovl_site.group_jpn[[3]])[names(ovl_site.group_jpn[[3]])=='Site']<-'site.group'
names(ovl_site.group_jpn[[2]])[names(ovl_site.group_jpn[[2]])=='Site']<-'site.group'
names(ovl_site.group_jpn[[1]])[names(ovl_site.group_jpn[[1]])=='Site']<-'site.group'

ovl_site.group_jpn[[3]]$site.group<-factor(ovl_site.group_jpn[[3]]$site.group,levels=levz)
ovl_site.group_jpn[[1]]$site.group<-factor(ovl_site.group_jpn[[1]]$site.group,levels=levz)
ovl_site.group_jpn[[2]]$site.group<-factor(ovl_site.group_jpn[[2]]$site.group,levels=levz)

#ovl_site.group_jpn[[1]]<-ovl_site.group_jpn[[1]][-which(is.na(ovl_site.group_jpn[[1]]$site.group)),]

ap1<-jpn_sp_site_pco[jpn_sp_site_pco$FG %in% c(15, 10, 8, 2,6, 12, 4,1, 16),]
ap1$FG<-factor(ap1$FG)

ap1$ThermalAffinity2<-factor(ap1$ThermalAffinity2, levels=c('tropical', 'subtropical'))
ovl_site.group_jpn[[3]]$ThermalAffinity2<-factor(ovl_site.group_jpn[[3]]$ThermalAffinity2, levels=c('tropical', 'subtropical'))

p4<-ppp+geom_point(data=ap1, aes(x=A1, y=A2, colour=ThermalAffinity2))+
  geom_sf(data=ovl_site.group_jpn[[3]], aes(fill=ThermalAffinity2), alpha=0.5)+
  geom_text(data=ovl_site.group_jpn[[1]], aes(x=0.9, y=0.9, label=paste(round(Freq, 2)*100,'%', sep='')))+
  geom_text(data=ovl_site.group_jpn[[2]], aes(x=0.7, y=0.7, label=paste(round(Freq, 2)*100,'%', sep='')), colour='red')+
  scale_y_continuous(paste('PC2', sprintf('(%0.1f%% explained var.)',
                                          100* func_dudi$eig[2]/sum(func_dudi$eig[func_dudi$eig>0.007]))))+
  scale_x_continuous(paste('PC1', sprintf('(%0.1f%% explained var.)',
                                          100* func_dudi$eig[1]/sum(func_dudi$eig[func_dudi$eig>0.007]))))+
  facet_grid(FG~site.group)+theme_minimal()+theme(legend.position = "none")

p4a<-ppp+
  geom_sf(data=ovl_site.group_jpn[[3]], aes(colour=ThermalAffinity2), fill=NA)+
  geom_text(data=ovl_site.group_jpn[[1]], aes(x=0.9, y=0.9, label=paste(round(Freq, 2)*100,'%', sep='')))+
  geom_text(data=ovl_site.group_jpn[[2]], aes(x=0.7, y=0.7, label=paste(round(Freq, 2)*100,'%', sep='')), colour='red')+
  scale_y_continuous(paste('PC2', sprintf('(%0.1f%% explained var.)',
                                          100* func_dudi$eig[2]/sum(func_dudi$eig[func_dudi$eig>0.007]))))+
  scale_x_continuous(paste('PC1', sprintf('(%0.1f%% explained var.)',
                                          100* func_dudi$eig[1]/sum(func_dudi$eig[func_dudi$eig>0.007]))))+
  facet_grid(FG~site.group)+theme_minimal()+theme(legend.position = "none")

### export Fig 4 componenets seperately for construction in CorelDraw ###

ggsave('C:/coral_fish/outputs/fig5_Australia_trends_pcoafix.eps',
       plot=p1,width = 21, height = 30, units = "cm")
ggsave('C:/coral_fish/outputs/fig5_Australia_niche_pcoafix.eps',
       plot=p2,width = 21, height = 30, units = "cm")
ggsave('C:/coral_fish/outputs/fig5_Japan_trends_pcoafix.eps',
       plot=p3,width = 21, height = 30, units = "cm")
ggsave('C:/coral_fish/outputs/fig5_Japan_niche_pcoafix.eps',
       plot=p4,width = 21, height = 30, units = "cm")

ggsave('C:/coral_fish/outputs/fig5_Australia_nicheALT_pcoafix.eps',
       plot=p2a,width = 21, height = 30, units = "cm")
ggsave('C:/coral_fish/outputs/fig5_Japan_nicheALT_pcoafix.eps',
       plot=p4a,width = 21, height = 30, units = "cm")
########################################################################

#### Calc biomass weighted centroid for each trop/temp niche ###
#### then calc distance between to complement functional overlap #####


# AUS by site
aus_site_cent<-aus_sp_site_pco%>%filter(FG %in% c(15, 10, 8, 2,6, 12, 4,1, 16))%>%group_by(Site, site.group, FG, ThermalAffinity2)%>%
  summarise(wtmn_A1=sum(A1*cor_biom)/sum(cor_biom), wtmn_A2=sum(A2*cor_biom)/sum(cor_biom), .groups='keep')

aus_site_dist<-NULL # euclidean dist function per SITE
for(j in unique(aus_site_cent$Site))
{
  s1<-aus_site_cent[aus_site_cent$Site==j,]
  
  for(k in unique(s1$FG))
  {
    s1_fg1<-s1[s1$FG==k,]
    if(nrow(s1_fg1)==1){out<-data.frame(s1_fg1[1,1:3], cent_d=NA)}else{
    out<-data.frame(s1_fg1[1,1:3], cent_d=as.numeric(dist(s1_fg1[,5:6])))}
    aus_site_dist<-rbind(aus_site_dist, out)
  }
}
aus_site_dist$FG<-as.character(aus_site_dist$FG)

# calc correlation
ovl_test_aus<-aus_trop_prop%>%group_by(FG, Site)%>%summarise(trop_met=mean(trop_met))%>%
  left_join(., ovl_site_aus[[1]][c(1,2,5)], by=c('FG', 'Site'))%>%
  left_join(., aus_site_dist[,c(1,3,4)], by=c('FG', 'Site'))

qplot(data=ovl_test_aus, x=Freq, y=cent_d)+facet_wrap(~FG, scales='free')
qplot(data=ovl_test_aus, x=cent_d, y=trop_met)+facet_wrap(~FG, scales='free')

# JPN by site
jpn_site_cent<-jpn_sp_site_pco%>%filter(FG %in% c(15, 10, 8, 2,6, 12, 4,1, 16))%>%group_by(Site, site.group, FG, ThermalAffinity2)%>%
  summarise(wtmn_A1=sum(A1*cor_biom)/sum(cor_biom), wtmn_A2=sum(A2*cor_biom)/sum(cor_biom), .groups='keep')

jpn_site_dist<-NULL # euclidean dist function per SITE
for(j in unique(jpn_site_cent$Site))
{
  s1<-jpn_site_cent[jpn_site_cent$Site==j,]
  
  for(k in unique(s1$FG))
  {
    s1_fg1<-s1[s1$FG==k,]
    if(nrow(s1_fg1)==1){out<-data.frame(s1_fg1[1,1:3], cent_d=NA)}else{
      out<-data.frame(s1_fg1[1,1:3], cent_d=as.numeric(dist(s1_fg1[,5:6])))}
    jpn_site_dist<-rbind(jpn_site_dist, out)
  }
}
jpn_site_dist$FG<-as.character(jpn_site_dist$FG)

# calc correlation
ovl_test_jpn<-jpn_trop_prop%>%group_by(FG, SiteID)%>%summarise(trop_met=mean(trop_met))%>%
  left_join(., ovl_site_jpn[[1]][c(1,2,5)], by=c('FG', 'SiteID'='Site'))%>%
  left_join(., jpn_site_dist[,c(1,3,4)], by=c('FG', 'SiteID'='Site'))

qplot(data=ovl_test_jpn, x=Freq, y=cent_d)+facet_wrap(~FG, scales='free')
qplot(data=ovl_test_jpn, x=cent_d, y=trop_met)+facet_wrap(~FG, scales='free')

# calc corr coef for both Aus and Jpn

# omit FG 16 as no competition, also 1(jpn),4(aus),6,12 due to too many NAs
kend_tests<-rbind(
  ovl_test_jpn%>%filter(trop_met!=0 & !FG%in%c(16, 1,6, 12))%>%group_by(FG)%>%summarise(comp_cent_p=cor.test(cent_d, Freq, method = 'kendall')$p.value,
                                                                         comp_cent_est=cor.test(cent_d, Freq, method = 'kendall')$estimate,
                                                                         cent_trop_p=cor.test(trop_met, cent_d, method = 'kendall')$p.value,
                                                                         cent_trop_est=cor.test(trop_met, cent_d, method = 'kendall')$estimate),
  ovl_test_aus%>%filter(trop_met!=0 & !FG%in%c(16, 4,6, 12))%>%group_by(FG)%>%summarise(comp_cent_p=cor.test(cent_d, Freq, method = 'kendall')$p.value,
                                                                         comp_cent_est=cor.test(cent_d, Freq, method = 'kendall')$estimate,
                                                                         cent_trop_p=cor.test(trop_met, cent_d, method = 'kendall')$p.value,
                                                                         cent_trop_est=cor.test(trop_met, cent_d, method = 'kendall')$estimate))

#write.csv(kend_tests, 'C:/tropicalisation_paper/reviewer_analyses/cent_overlap_correlation_tropicalization_pcoafix.csv', quote=F, row.names=F)

# now write out summary vals
summr<-rbind(
  ovl_test_jpn%>%filter(trop_met!=0)%>%group_by(FG)%>%summarise(comp_cent_mean=mean(cent_d, na.rm=T), comp_cent_sd=sd(cent_d, na.rm=T)),
  ovl_test_aus%>%filter(trop_met!=0)%>%group_by(FG)%>%summarise(comp_cent_mean=mean(cent_d, na.rm=T), comp_cent_sd=sd(cent_d, na.rm=T)))
                                                                                                                    
#write.csv(summr, 'C:/tropicalisation_paper/reviewer_analyses/cent_meanz_summary_pcoafix.csv', quote=F, row.names=F)

# make lat plots
## AUSTRALIA
aus_ovl_comp<-left_join(ovl_site_aus[[1]],aus_sp_site%>%group_by(Site)%>%
                          summarise_all(first)%>%dplyr::select(Site, Lat), by='Site')

aus_site_dist_lat<-left_join(aus_site_dist, aus_ovl_comp[c('Site', 'Lat')], by='Site')


#aus_ovl_comp<-aus_ovl_comp[-which(is.na(aus_ovl_comp$Site)),] # remove na row from FG 16 no overlap

aus_trop_prop$FG<-factor(aus_trop_prop$FG,
                         levels=c(15, 10, 8, 2,6, 12, 4,1, 16))
aus_ovl_comp$FG<-factor(aus_ovl_comp$FG,
                         levels=c(15, 10, 8, 2,6, 12, 4,1, 16))
aus_site_dist_lat$FG<-factor(aus_site_dist_lat$FG,
                        levels=c(15, 10, 8, 2,6, 12, 4,1, 16))

sg_lat_spans<-data.frame(xmin=c(-23.4, -24.8, -26.61, -28.19, -29.9, -29.97), 
                         xmax=c(-24.1, -25.3, -26.98, -28.616, -30.96, -30.3))
p_aus_cent<-ggplot()+
  geom_rect(data=sg_lat_spans, aes(ymin=0, ymax=2, xmin=xmin, xmax=xmax), fill='grey', alpha=0.5)+
  geom_smooth(data=aus_trop_prop, aes(x=Lat, y=trop_met), se=F, colour='black')+
  geom_smooth(data=aus_ovl_comp, aes(x=Lat, y=Freq), colour='red', se=F)+
  geom_smooth(data=aus_site_dist_lat%>%filter(!FG%in%c(16, 4,6, 12)), aes(x=Lat, y=cent_d), colour='purple', se=F)+
  geom_point(data=filter(aus_trop_prop, trop_met<2), aes(x=Lat, y=trop_met))+
  geom_point(data=aus_ovl_comp, aes(x=Lat, y=Freq), colour='red')+
  geom_point(data=aus_site_dist_lat, aes(x=Lat, y=cent_d), colour='purple')+
  theme_bw()+facet_grid(FG~.)+scale_x_reverse()

## JAPAN
jpn_ovl_comp<-left_join(ovl_site_jpn[[1]],jpn_sp_site%>%group_by(SiteID)%>%
                          summarise_all(first)%>%dplyr::select(SiteID, lat), by=c('Site'='SiteID'))

jpn_site_dist_lat<-left_join(jpn_site_dist, jpn_ovl_comp[c('Site', 'lat')], by='Site')


#jpn_ovl_comp<-jpn_ovl_comp[-which(is.na(jpn_ovl_comp$Site)),] # remove na row from FG 16 no overlap

jpn_trop_prop$FG<-factor(jpn_trop_prop$FG,
                         levels=c(15, 10, 8, 2,6, 12, 4,1, 16))
jpn_ovl_comp$FG<-factor(jpn_ovl_comp$FG,
                        levels=c(15, 10, 8, 2,6, 12, 4,1, 16))
jpn_site_dist_lat$FG<-factor(jpn_site_dist_lat$FG,
                             levels=c(15, 10, 8, 2,6, 12, 4,1, 16))

sg_lat_spans<-data.frame(xmin=c(24.2, 26.2, 28.5, 31,   32.7,  33.38, 34.6), 
                         xmax=c(24.5, 28.4, 30.5, 31.6, 32.82, 33.5, 35 ))

p_jpn_cent<-ggplot()+
  geom_rect(data=sg_lat_spans, aes(ymin=0, ymax=2, xmin=xmin, xmax=xmax), fill='grey', alpha=0.5)+
  geom_smooth(data=jpn_trop_prop, aes(x=lat, y=trop_met), se=F, colour='black')+
  geom_smooth(data=jpn_ovl_comp, aes(x=lat, y=Freq), colour='red', se=F)+
  geom_smooth(data=jpn_site_dist_lat%>%filter(!FG%in%c(16, 1,6, 12)), aes(x=lat, y=cent_d), colour='purple', se=F)+
  geom_point(data=filter(jpn_trop_prop, trop_met<2), aes(x=lat, y=trop_met))+
  geom_point(data=jpn_ovl_comp, aes(x=lat, y=Freq), colour='red')+
  geom_point(data=jpn_site_dist_lat, aes(x=lat, y=cent_d), colour='purple')+
  theme_bw()+facet_grid(FG~.)

#check against plot (site.group)

#AUSTRALIA
aus_site.group_cent<-aus_sp_site_pco%>%filter(FG %in% c(15, 10, 8, 2,6, 12, 4,1, 16))%>%group_by(site.group, FG, ThermalAffinity2)%>%
  summarise(wtmn_A1=sum(A1*cor_biom)/sum(cor_biom), wtmn_A2=sum(A2*cor_biom)/sum(cor_biom), .groups='keep')

aus_site.group_dist<-NULL # euclidean dist function per SITE.GROUP
for(j in unique(aus_site.group_cent$site.group))
{
  s1<-aus_site.group_cent[aus_site.group_cent$site.group==j,]
  
  for(k in unique(s1$FG))
  {
    s1_fg1<-s1[s1$FG==k,]
    if(nrow(s1_fg1)==1){out<-data.frame(s1_fg1[1,1:2], cent_d=NA)}else{
      out<-data.frame(s1_fg1[1,1:2], cent_d=as.numeric(dist(s1_fg1[,4:5])))}
    aus_site.group_dist<-rbind(aus_site.group_dist, out)
  }
}

levz2<-c('trop.base', 'trans.bay', 'trans.offshore', 'trans.temp',
         'temp.offshore','temp.inshore')

aus_site.group_dist$site.group<-factor(aus_site.group_dist$site.group,levels=levz2)
aus_site.group_dist$FG<-factor(aus_site.group_dist$FG, levels = c(15, 10, 8, 2, 6, 12, 4, 1, 16))

aus_site.group_cent$FG<-factor(aus_site.group_cent$FG, levels = c(15, 10, 8, 2, 6, 12, 4, 1, 16))
aus_site.group_cent$site.group<-factor(aus_site.group_dist$site.group,levels=levz2)

ovl_site.group_aus[[3]]$FG<-factor(ovl_site.group_aus[[3]]$FG,levels=c(15, 10, 8, 2, 6, 12, 4, 1, 16))
ovl_site.group_aus[[1]]$FG<-factor(ovl_site.group_aus[[1]]$FG,levels=c(15, 10, 8, 2, 6, 12, 4, 1, 16))

ovl_site.group_aus[[3]]$site.group<-factor(ovl_site.group_aus[[3]]$site.group,levels=levz2)
ovl_site.group_aus[[1]]$site.group<-factor(ovl_site.group_aus[[1]]$site.group,levels=levz2)

aus_trait_cent<-ppp+geom_sf(data=ovl_site.group_aus[[3]], aes(colour=ThermalAffinity2), fill=NA)+
  scale_y_continuous(paste('PC2', sprintf('(%0.1f%% explained var.)',
                                          100* func_dudi$eig[2]/sum(func_dudi$eig[func_dudi$eig>0.007]))))+
  scale_x_continuous(paste('PC1', sprintf('(%0.1f%% explained var.)',
                                          100* func_dudi$eig[1]/sum(func_dudi$eig[func_dudi$eig>0.007]))))+
  facet_grid(FG~site.group)+theme_minimal()+theme(legend.position = "none")+
  geom_text(data=ovl_site.group_aus[[1]], aes(x=0.7, y=0.7, label=paste(round(Freq, 2)*100,'%', sep='')), colour='red')+
  geom_text(data=aus_site.group_dist, aes(x=-0.7, y=-0.7, label= round((aus_site.group_dist$cent_d*10), 1)), colour='purple')+
  geom_point(data=aus_site.group_cent, aes(x=wtmn_A1 , y=wtmn_A2, colour=ThermalAffinity2))

# JAPAN

#check against plot (site.group)
jpn_site.group_cent<-jpn_sp_site_pco%>%filter(FG %in% c(15, 10, 8, 2,6, 12, 4,1, 16))%>%group_by(site.group, FG, ThermalAffinity2)%>%
  summarise(wtmn_A1=sum(A1*cor_biom)/sum(cor_biom), wtmn_A2=sum(A2*cor_biom)/sum(cor_biom), .groups='keep')

jpn_site.group_dist<-NULL # euclidean dist function per SITE.GROUP
for(j in unique(jpn_site.group_cent$site.group))
{
  s1<-jpn_site.group_cent[jpn_site.group_cent$site.group==j,]
  
  for(k in unique(s1$FG))
  {
    s1_fg1<-s1[s1$FG==k,]
    if(nrow(s1_fg1)==1){out<-data.frame(s1_fg1[1,1:2], cent_d=NA)}else{
      out<-data.frame(s1_fg1[1,1:2], cent_d=as.numeric(dist(s1_fg1[,4:5])))}
    jpn_site.group_dist<-rbind(jpn_site.group_dist, out)
  }
}

levz2<-c('trop.base', 'trop.island', 'trans.island', 
        'trans.headld', 'trans.inland','temp.headld')

jpn_site.group_dist$site.group<-factor(jpn_site.group_dist$site.group,levels=levz2)
jpn_site.group_dist$FG<-factor(jpn_site.group_dist$FG, levels = c(15, 10, 8, 2, 6, 12, 4, 1, 16))

jpn_site.group_cent$FG<-factor(jpn_site.group_cent$FG, levels = c(15, 10, 8, 2, 6, 12, 4, 1, 16))
jpn_site.group_cent$site.group<-factor(jpn_site.group_dist$site.group,levels=levz2)

ovl_site.group_jpn[[3]]$FG<-factor(ovl_site.group_jpn[[3]]$FG,levels=c(15, 10, 8, 2, 6, 12, 4, 1, 16))
ovl_site.group_jpn[[1]]$FG<-factor(ovl_site.group_jpn[[1]]$FG,levels=c(15, 10, 8, 2, 6, 12, 4, 1, 16))

ovl_site.group_jpn[[3]]$site.group<-factor(ovl_site.group_jpn[[3]]$site.group,levels=levz2)
ovl_site.group_jpn[[1]]$site.group<-factor(ovl_site.group_jpn[[1]]$site.group,levels=levz2)

jpn_trait_cent<-ppp+geom_sf(data=ovl_site.group_jpn[[3]], aes(colour=ThermalAffinity2), fill=NA)+
  scale_y_continuous(paste('PC2', sprintf('(%0.1f%% explained var.)',
                                          100* func_dudi$eig[2]/sum(func_dudi$eig[func_dudi$eig>0.007]))))+
  scale_x_continuous(paste('PC1', sprintf('(%0.1f%% explained var.)',
                                          100* func_dudi$eig[1]/sum(func_dudi$eig[func_dudi$eig>0.007]))))+
  facet_grid(FG~site.group)+theme_minimal()+theme(legend.position = "none")+
  geom_text(data=ovl_site.group_jpn[[1]], aes(x=0.7, y=0.7, label=paste(round(Freq, 2)*100,'%', sep='')), colour='red')+
  geom_text(data=jpn_site.group_dist, aes(x=-0.7, y=-0.7, label= round((jpn_site.group_dist$cent_d*10), 1)), colour='purple')+
  geom_point(data=jpn_site.group_cent, aes(x=wtmn_A1 , y=wtmn_A2, colour=ThermalAffinity2))

### export Fig 4 componenets seperately for construction in CorelDraw ###

ggsave('C:/tropicalisation_paper/reviewer_analyses/fig_Australia_trends_cent.eps',
       plot=p_aus_cent,width = 21, height = 30, units = "cm")

ggsave('C:/tropicalisation_paper/reviewer_analyses/fig_Japan_trends_cent.eps',
       plot=p_jpn_cent,width = 21, height = 30, units = "cm")

ggsave('C:/tropicalisation_paper/reviewer_analyses/fig_Australia_nicheALT_cent.eps',
       plot=aus_trait_cent,width = 21, height = 30, units = "cm")
ggsave('C:/tropicalisation_paper/reviewer_analyses/fig_Japan_nicheALT_cent.eps',
       plot=jpn_trait_cent,width = 21, height = 30, units = "cm")

########################################################################

#### Calculate overlap at site level over first four ###
#### PCoA axis and combine based on weighted average (variance explained by each PCoA pair) #####

# run function per site for lat plots with tropicalization metric
ovl_site_aus1<-funcOvl(pcodat=aus_sp_site_pco, FGz=c(15, 10, 8, 2,6, 12, 4,1, 16), bywhat='site', mkern=F, selPC1="A1", selPC2="A2")
ovl_site_aus2<-funcOvl(pcodat=aus_sp_site_pco, FGz=c(15, 10, 8, 2,6, 12, 4,1, 16), bywhat='site', mkern=F, selPC1="A1", selPC2="A3")
ovl_site_aus3<-funcOvl(pcodat=aus_sp_site_pco, FGz=c(15, 10, 8, 2,6, 12, 4,1, 16), bywhat='site', mkern=F, selPC1="A1", selPC2="A4")
ovl_site_aus4<-funcOvl(pcodat=aus_sp_site_pco, FGz=c(15, 10, 8, 2,6, 12, 4,1, 16), bywhat='site', mkern=F, selPC1="A2", selPC2="A3")
ovl_site_aus5<-funcOvl(pcodat=aus_sp_site_pco, FGz=c(15, 10, 8, 2,6, 12, 4,1, 16), bywhat='site', mkern=F, selPC1="A2", selPC2="A4")
ovl_site_aus6<-funcOvl(pcodat=aus_sp_site_pco, FGz=c(15, 10, 8, 2,6, 12, 4,1, 16), bywhat='site', mkern=F, selPC1="A3", selPC2="A4")

ovl_site_jpn1<-funcOvl(pcodat=jpn_sp_site_pco, FGz=c(15, 10, 8, 2,6, 12, 4,1, 16), bywhat='site', mkern=F, selPC1="A1", selPC2="A2")
ovl_site_jpn2<-funcOvl(pcodat=jpn_sp_site_pco, FGz=c(15, 10, 8, 2,6, 12, 4,1, 16), bywhat='site', mkern=F, selPC1="A1", selPC2="A3")
ovl_site_jpn3<-funcOvl(pcodat=jpn_sp_site_pco, FGz=c(15, 10, 8, 2,6, 12, 4,1, 16), bywhat='site', mkern=F, selPC1="A1", selPC2="A4")
ovl_site_jpn4<-funcOvl(pcodat=jpn_sp_site_pco, FGz=c(15, 10, 8, 2,6, 12, 4,1, 16), bywhat='site', mkern=F, selPC1="A2", selPC2="A3")
ovl_site_jpn5<-funcOvl(pcodat=jpn_sp_site_pco, FGz=c(15, 10, 8, 2,6, 12, 4,1, 16), bywhat='site', mkern=F, selPC1="A2", selPC2="A4")
ovl_site_jpn6<-funcOvl(pcodat=jpn_sp_site_pco, FGz=c(15, 10, 8, 2,6, 12, 4,1, 16), bywhat='site', mkern=F, selPC1="A3", selPC2="A4")

# fill overlap with 0's where no overlap
ovl_site_aus1[[1]]<-ovl_site_aus1[[1]]%>%tidyr::complete(Site, FG, fill=list(Freq=0))
ovl_site_aus1[[2]]<-ovl_site_aus1[[2]]%>%tidyr::complete(Site, FG, fill=list(Freq=0))
ovl_site_aus2[[1]]<-ovl_site_aus2[[1]]%>%tidyr::complete(Site, FG, fill=list(Freq=0))
ovl_site_aus2[[2]]<-ovl_site_aus2[[2]]%>%tidyr::complete(Site, FG, fill=list(Freq=0))
ovl_site_aus3[[1]]<-ovl_site_aus3[[1]]%>%tidyr::complete(Site, FG, fill=list(Freq=0))
ovl_site_aus3[[2]]<-ovl_site_aus3[[2]]%>%tidyr::complete(Site, FG, fill=list(Freq=0))
ovl_site_aus4[[1]]<-ovl_site_aus4[[1]]%>%tidyr::complete(Site, FG, fill=list(Freq=0))
ovl_site_aus4[[2]]<-ovl_site_aus4[[2]]%>%tidyr::complete(Site, FG, fill=list(Freq=0))
ovl_site_aus5[[1]]<-ovl_site_aus5[[1]]%>%tidyr::complete(Site, FG, fill=list(Freq=0))
ovl_site_aus5[[2]]<-ovl_site_aus5[[2]]%>%tidyr::complete(Site, FG, fill=list(Freq=0))
ovl_site_aus6[[1]]<-ovl_site_aus6[[1]]%>%tidyr::complete(Site, FG, fill=list(Freq=0))
ovl_site_aus6[[2]]<-ovl_site_aus6[[2]]%>%tidyr::complete(Site, FG, fill=list(Freq=0))

ovl_site_jpn1[[1]]<-ovl_site_jpn1[[1]]%>%tidyr::complete(Site, FG, fill=list(Freq=0))
ovl_site_jpn1[[2]]<-ovl_site_jpn1[[2]]%>%tidyr::complete(Site, FG, fill=list(Freq=0))
ovl_site_jpn2[[1]]<-ovl_site_jpn2[[1]]%>%tidyr::complete(Site, FG, fill=list(Freq=0))
ovl_site_jpn2[[2]]<-ovl_site_jpn2[[2]]%>%tidyr::complete(Site, FG, fill=list(Freq=0))
ovl_site_jpn3[[1]]<-ovl_site_jpn3[[1]]%>%tidyr::complete(Site, FG, fill=list(Freq=0))
ovl_site_jpn3[[2]]<-ovl_site_jpn3[[2]]%>%tidyr::complete(Site, FG, fill=list(Freq=0))
ovl_site_jpn4[[1]]<-ovl_site_jpn4[[1]]%>%tidyr::complete(Site, FG, fill=list(Freq=0))
ovl_site_jpn4[[2]]<-ovl_site_jpn4[[2]]%>%tidyr::complete(Site, FG, fill=list(Freq=0))
ovl_site_jpn5[[1]]<-ovl_site_jpn5[[1]]%>%tidyr::complete(Site, FG, fill=list(Freq=0))
ovl_site_jpn5[[2]]<-ovl_site_jpn5[[2]]%>%tidyr::complete(Site, FG, fill=list(Freq=0))
ovl_site_jpn6[[1]]<-ovl_site_jpn6[[1]]%>%tidyr::complete(Site, FG, fill=list(Freq=0))
ovl_site_jpn6[[2]]<-ovl_site_jpn6[[2]]%>%tidyr::complete(Site, FG, fill=list(Freq=0))

# combine all metrics together
aus_all_comp<-data.frame(ovl_site_aus1[[1]][,c(1,2,5)], ovl_site_aus2[[1]][,5],  ovl_site_aus3[[1]][,5],
                         ovl_site_aus4[[1]][,5],  ovl_site_aus5[[1]][,5],  ovl_site_aus6[[1]][,5])
                      
aus_all_filt<-data.frame(ovl_site_aus1[[2]][,c(1,2,5)], ovl_site_aus2[[2]][,5],  ovl_site_aus3[[2]][,5],
                         ovl_site_aus4[[2]][,5],  ovl_site_aus5[[2]][,5],  ovl_site_aus6[[2]][,5])

jpn_all_comp<-data.frame(ovl_site_jpn1[[1]][,c(1,2,5)], ovl_site_jpn2[[1]][,5],  ovl_site_jpn3[[1]][,5],
                         ovl_site_jpn4[[1]][,5],  ovl_site_jpn5[[1]][,5],  ovl_site_jpn6[[1]][,5])

jpn_all_filt<-data.frame(ovl_site_jpn1[[2]][,c(1,2,5)], ovl_site_jpn2[[2]][,5],  ovl_site_jpn3[[2]][,5],
                         ovl_site_jpn4[[2]][,5],  ovl_site_jpn5[[2]][,5],  ovl_site_jpn6[[2]][,5])
                      
# plot to visualise how similar/different metrics calcualted from different PCoA combos are
aus_ovl_comp<-left_join(aus_all_comp,aus_sp_site%>%group_by(Site)%>%
                          summarise_all(first)%>%dplyr::select(Site, Lat), by='Site')
aus_ovl_filt<-left_join(aus_all_filt,aus_sp_site%>%group_by(Site)%>%
                          summarise_all(first)%>%dplyr::select(Site, Lat), by='Site')

aus_ovl_comp<-pivot_longer(aus_ovl_comp, cols=Freq:Freq.5)
aus_ovl_filt<-pivot_longer(aus_ovl_filt, cols=Freq:Freq.5)

sg_lat_spans<-data.frame(xmin=c(-23.4, -24.8, -26.61, -28.19, -29.9, -29.97), 
                         xmax=c(-24.1, -25.3, -26.98, -28.616, -30.96, -30.3))
ggplot()+
  geom_rect(data=sg_lat_spans, aes(ymin=0, ymax=2, xmin=xmin, xmax=xmax), fill='grey', alpha=0.5)+
  geom_smooth(data=aus_trop_prop, aes(x=Lat, y=trop_met), se=F, colour='black')+
  geom_smooth(data=aus_ovl_comp, aes(x=Lat, y=value, group=name), colour='red', se=F)+
  geom_smooth(data=aus_ovl_filt, aes(x=Lat, y=value, group=name), colour='blue', se=F)+
  geom_point(data=filter(aus_trop_prop, trop_met<2), aes(x=Lat, y=trop_met))+
  geom_point(data=aus_ovl_comp, aes(x=Lat, y=value), colour='red')+
  geom_point(data=aus_ovl_filt, aes(x=Lat, y=value), colour='blue')+
  theme_bw()+facet_grid(FG~.)+scale_x_reverse()

#var explained by PCs 1:4
my_var_exp<-100* func_dudi$eig[1:4]/sum(func_dudi$eig[func_dudi$eig>0.007])
# get combinations 
f1<-my_var_exp[1]+my_var_exp[2] # selPC1="A1", selPC2="A2"
f2<-my_var_exp[1]+my_var_exp[3] # selPC1="A1", selPC2="A3"
f3<-my_var_exp[1]+my_var_exp[4] # selPC1="A1", selPC2="A4"
f4<-my_var_exp[2]+my_var_exp[3] # selPC1="A2", selPC2="A3"
f5<-my_var_exp[2]+my_var_exp[4] # selPC1="A2", selPC2="A4"
f6<-my_var_exp[3]+my_var_exp[4] # selPC1="A3", selPC2="A4"

# Reviewer suggestion to identify any axes (dimensions where we have zero overlap
# in which case there is no overlap at all and average shoudl instead yield 0
fun2<-function(x){which(as.numeric(x[c("Freq", "Freq.1", "Freq.2", "Freq.3", "Freq.4", "Freq.5")])==0)}

aus_all_comp$n_zero<-apply(aus_all_comp,1, FUN=fun2)
aus_all_filt$n_zero<-apply(aus_all_filt,1, FUN=fun2)
jpn_all_comp$n_zero<-apply(jpn_all_comp,1, FUN=fun2)
jpn_all_filt$n_zero<-apply(jpn_all_filt,1, FUN=fun2)

# wherever there is a zero, replace with NA. This means that in weighted.mean we get NA overlap if
# any NA occurs in 6 pairs of overlap
aus_all_comp[aus_all_comp==0]<-NA
aus_all_filt[aus_all_filt==0]<-NA
jpn_all_comp[jpn_all_comp==0]<-NA
jpn_all_filt[jpn_all_filt==0]<-NA

fun1<-function(x){weighted.mean(as.numeric(x[c("Freq", "Freq.1", "Freq.2", "Freq.3", "Freq.4", "Freq.5")]), c(f1, f2, f3, f4, f5, f6))}

aus_all_comp$Freq_WM<-apply(aus_all_comp,1, FUN=fun1)
aus_all_filt$Freq_WM<-apply(aus_all_filt,1, FUN=fun1)

jpn_all_comp$Freq_WM<-apply(jpn_all_comp,1, FUN=fun1)
jpn_all_filt$Freq_WM<-apply(jpn_all_filt,1, FUN=fun1)

#change NA Freq_WM back to zero
aus_all_comp[is.na(aus_all_comp$Freq_WM),]$Freq_WM<-0
aus_all_filt[is.na(aus_all_filt$Freq_WM),]$Freq_WM<-0
jpn_all_comp[is.na(jpn_all_comp$Freq_WM),]$Freq_WM<-0
jpn_all_filt[is.na(jpn_all_filt$Freq_WM),]$Freq_WM<-0

#### Calc correlation coef & sig between functional ###
#### niche metrics and tropicalisation footprints #####

ovl_test_jpn<-jpn_trop_prop%>%group_by(FG, SiteID)%>%summarise(trop_met=mean(trop_met))%>%
  left_join(., jpn_all_comp[c("Site", "FG", "Freq_WM")], by=c('FG', 'SiteID'='Site'))%>%
  left_join(., jpn_all_filt[c("Site", "FG", "Freq_WM")], by=c('FG', 'SiteID'='Site'))

ovl_test_aus<-aus_trop_prop%>%group_by(FG, Site)%>%summarise(trop_met=mean(trop_met))%>%
  left_join(., aus_all_comp[c("Site", "FG", "Freq_WM")], by=c('FG', 'Site'))%>%
  left_join(., aus_all_filt[c("Site", "FG", "Freq_WM")], by=c('FG', 'Site'))

qplot(data=ovl_test_jpn, x=Freq_WM.y, y=trop_met)+facet_wrap(~FG, scales='free')
qplot(data=ovl_test_aus, x=Freq_WM.y, y=trop_met)+facet_wrap(~FG, scales='free')

# tests for aus and japan then write out
# omit FG 16 as no competition
kend_tests<-rbind(
  ovl_test_jpn%>%filter(trop_met!=0 & FG!=16)%>%group_by(FG)%>%summarise(comp_p=cor.test(trop_met, Freq_WM.x, method = 'kendall')$p.value,
                                                                comp_est=cor.test(trop_met, Freq_WM.x, method = 'kendall')$estimate,
                                                                filt_p=cor.test(trop_met, Freq_WM.y, method = 'kendall')$p.value,
                                                                filt_est=cor.test(trop_met, Freq_WM.y, method = 'kendall')$estimate),
  ovl_test_aus%>%filter(trop_met!=0 & FG!=16)%>%group_by(FG)%>%summarise(comp_p=cor.test(trop_met, Freq_WM.x, method = 'kendall')$p.value,
                                                                comp_est=cor.test(trop_met, Freq_WM.x, method = 'kendall')$estimate,
                                                                filt_p=cor.test(trop_met, Freq_WM.y, method = 'kendall')$p.value,
                                                                filt_est=cor.test(trop_met, Freq_WM.y, method = 'kendall')$estimate))
#write.csv(kend_tests, 'C:/coral_fish/outputs/func_overlap_correlation_tropicalization_pcoafix.csv', quote=F, row.names=F)

# now write out summary vals
summr<-rbind(
  ovl_test_jpn%>%filter(trop_met!=0)%>%group_by(FG)%>%summarise(comp_mean=mean(Freq_WM.x, na.rm=T), comp_sd=sd(Freq_WM.x, na.rm=T),
                                                                filt_mean=mean(Freq_WM.y), filt_sd=sd(Freq_WM.y)),
  ovl_test_aus%>%filter(trop_met!=0)%>%group_by(FG)%>%summarise(comp_mean=mean(Freq_WM.x), comp_sd=sd(Freq_WM.x),
                                                                filt_mean=mean(Freq_WM.y), filt_sd=sd(Freq_WM.y)))
#write.csv(summr, 'C:/coral_fish/outputs/func_overlap_meanz_summary_pcoafix.csv', quote=F, row.names=F)


# Do correlation test between competition and subtropical biomass

ovl_test_jpn<-left_join(ovl_test_jpn, bio_jpn%>%filter(.,ThermalAffinity2=='subtropical')%>%
                          group_by(FG, SiteID)%>%summarise(cor_biom=mean(cor_biom)), by=c('FG', 'SiteID'))

ovl_test_aus<-left_join(ovl_test_aus, bio_aus%>%filter(.,ThermalAffinity2=='subtropical')%>%
                          group_by(FG, Site)%>%summarise(cor_biom=mean(cor_biom)), by=c('FG', 'Site'))

qplot(data=ovl_test_jpn, x=Freq_WM.x, y=cor_biom^0.25)+facet_wrap(~FG, scales='free')
qplot(data=ovl_test_aus, x=Freq_WM.x, y=cor_biom^0.25)+facet_wrap(~FG, scales='free')

ovl_test_aus%>%filter(trop_met!=0& FG!=16)%>%group_by(FG)%>%summarise(comp_p=cor.test(Freq_WM.x, cor_biom^0.25, method = 'kendall')$p.value,
                                                              comp_est=cor.test(Freq_WM.x, cor_biom^0.25, method = 'kendall')$estimate)

ovl_test_jpn%>%filter(trop_met!=0& FG!=16)%>%group_by(FG)%>%summarise(comp_p=cor.test(Freq_WM.x, cor_biom^0.25, method = 'kendall')$p.value,
                                                              comp_est=cor.test(Freq_WM.x, cor_biom^0.25, method = 'kendall')$estimate)
# subtropical biomass sig positively correlated (p<0.05) with
# functional niche overlap in all FGs in both regions. apart from FG 6 in AUS and FG1 in JPN

## A tibble: 8 x 3
#FG       comp_p comp_est
#<chr>     <dbl>    <dbl>
#  1 1     0.00397      0.458
#2 10    0.0000308    0.567
#3 12    0.0000220    1    
#4 15    0.000461     0.558
#5 2     0.00721      0.426
#6 4     0.00149      0.730
#7 6     0.140        0.327
#8 8     0.000694     0.502

# A tibble: 8 x 3
#FG         comp_p comp_est
#<chr>       <dbl>    <dbl>
#  1 1     0.144          0.298
#2 10    0.000000708    0.728
#3 12    0.00000459     1    
#4 15    0.000000242    0.696
#5 2     0.00000558     0.687
#6 4     0.0108         0.405
#7 6     0.000000547    0.962
#8 8     0.0000145      0.652

##########################################################################################################################

#################### Paper Table 1 summary analyses ##########################
##############################################################################

# read in table 1 for summary analysis, populated with outputs from code in excel 
table1<-read.csv('table1_analysis.csv')
#join with thermal mean per FG data
tmean<-read.csv('thermal_midpoint_means.csv')
table1<-left_join(table1, tmean[,1:3], by=c('FG'='groupk19', 'region'))

table1$latclass2<-table1$latclass
table1[table1$latclass2=='above',]$latclass2<-'expected'
table1$latclass2<-  factor(table1$latclass2,
                           levels=c('well.below', 'below', 'expected'), ordered=T)

table1$Bdiff<-table1$max_temperateB/table1$max_tropicalB
table1$Bdiff2<-(table1$max_temperateB^0.25)/(table1$max_tropicalB^0.25)
qplot(data=table1, x=Bdiff, y=Bdiff2)

table(table1$latclass2, table1$diet)
table(table1$latclass2, table1$position)
table(table1$latclass2, table1$region)

#check significance of categorical variables
# all simulate p vals above exact so confident all are NS
chisq.test(table1$latclass2, table1$region, simulate.p.value = F)
#X-squared = 1.4, df = 2, p-value = 0.4966
chisq.test(table1$latclass2, table1$diet, simulate.p.value = F)
#X-squared = 10.05, df = 10, p-value = 0.4361
chisq.test(table1$latclass2, table1$position, simulate.p.value = F)
#X-squared = 4.3875, df = 6, p-value = 0.6244

# mod tropical biomass
m_tropB<-lm(max_tropicalB_4rt~latclass2, data=table1)
m_tropB<-glm(sqrt(max_tropicalB)~latclass2, data=table1, family=Gamma(link = "log"))
plot(m_tropB)
anova(m_tropB, test="F") #latclass2  2 0.082488        15     8.9162 0.0746 0.9284
emmeans(m_tropB, 'latclass2', type='response')

# mod temperate biomass
m_tempB_inCor<-lm(max_temperateB_4rt~latclass2, data=table1)
m_tempB<-glm(sqrt(max_temperateB)~latclass2, data=table1[table1$diet!='corallivores',],family=Gamma(link = "log"))
plot(m_tempB)
anova(m_tempB_inCor) #latclass2  2 14.047  7.0234  7.0534 0.00693 ** with corallivores
anova(m_tempB, test="F")# latclass2  2   3.2969        13     9.6903  0.02863 *
m_tempB_em<-emmeans(m_tempB, 'latclass2', type='response')

# mod Bdiff2 
#m_bdiff<-lm(Bdiff2~latclass2, data=table1)
m_tempB<-glm(sqrt(Bdiff)~latclass2, data=table1[table1$diet!='corallivores',],family=Gamma(link = "log"))
plot(m_bdiff)
anova(m_bdiff) #latclass2  2 2.3325 1.16622  9.7815 0.00191 **
emmeans(m_bdiff, 'latclass2')

# mod func_niceO 
m_funcO<-glm(cbind(func_niceO, 100-func_niceO)~latclass2, data=table1[table1$diet!='corallivores',], family='quasibinomial')
plot(m_funcO)
anova(m_funcO, test="F")# latclass2  2   157.99        13     128.34 8.5406 0.004283 **
m_funcO_em<-emmeans(m_funcO, 'latclass2', type='response')

# mod trop_funcA 
m_funcA<-glm(cbind(trop_funcA, 100-trop_funcA)~latclass2, data=table1, family='quasibinomial')
plot(m_funcA)
anova(m_funcA, test="F") #latclass2  2  0.80784        15     36.739 0.1643   0.85

# mod nspecies
m_nsp<-glm(nspecies~latclass2, data=table1, family='quasipoisson')
plot(m_nsp)
anova(m_nsp, test="F")# latclass2  2   142.45        15     180.99 6.0089 0.01211 *
m_nsp_em<-emmeans(m_nsp, 'latclass2', type='response')
sum((resid(m_nsp, type="pearson")^2))/df.residual(m_nsp)

# mod thermal mean 
m_therm<-lm(thermal_mn~latclass2, data=table1)
plot(m_therm)
anova(m_therm)# latclass2  2 0.020399 0.0101995  1.1382 0.3466
emmeans(m_therm, 'latclass2')

#png('C:/coral_fish/outputs/poleward_advance_correlation.png',width =8, height =8 , units ="in", res =300)
ggpairs(table1[c('max_tropicalB', 'max_temperateB','nspecies',
                 'func_niceO', 'trop_funcA', 'thermal_mn', 'latclass2')])
#dev.off()

# make plots for supplementary Fig S5 #

p1<-ggplot(data=table1, aes(x=latclass2, y=max_tropicalB))+
  geom_jitter(aes(colour=region),shape=1, size=2, width=0.1)+theme_bw()+xlab('FG poleward advance')+ ylab('Tropical biomass')+theme(legend.position='none') 
p2<-ggplot(data=table1, aes(x=latclass2, y=max_temperateB))+
  geom_jitter(aes(colour=region),shape=1, size=2, width=0.1)+theme_bw()+xlab('FG poleward advance')+ ylab('Temperate biomass')+theme(legend.position='none')+
  geom_pointrange(data=as.data.frame(m_tempB_em), aes(y=response, ymin=asymp.LCL, ymax=asymp.UCL), alpha=0.9)

ggplot(data=table1, aes(x=latclass2, y=Bdiff, colour=region))+geom_point(shape=1)+theme_bw()

p3<-ggplot(data=table1, aes(x=latclass2, y=func_niceO/100))+
  geom_jitter(aes(colour=region),shape=1, size=2, width=0.1)+theme_bw()+xlab('FG poleward advance')+ ylab('Tropical-temperate functional overlap')+
  theme(legend.position='none')+geom_pointrange(data=as.data.frame(m_funcO_em), aes(y=prob, ymin=asymp.LCL, ymax=asymp.UCL), alpha=0.9)

p4<-ggplot(data=table1, aes(x=latclass2, y=trop_funcA))+
  geom_jitter(aes(colour=region),shape=1, size=2, width=0.1)+theme_bw()+xlab('FG poleward advance')+ ylab('Tropical functional area maintained')+
  theme(legend.position='none')

p5<-ggplot(data=table1, aes(x=latclass2, y=nspecies))+
  geom_jitter(aes(colour=region),shape=1, size=2, width=0.1)+theme_bw()+xlab('FG poleward advance')+ ylab('Species richness')+
  theme(legend.position='none')+geom_pointrange(data=as.data.frame(m_nsp_em), aes(y=rate, ymin=asymp.LCL, ymax=asymp.UCL), alpha=0.9)

p6<-ggplot(data=table1, aes(x=latclass2, y=thermal_mn))+
  geom_jitter(aes(colour=region),shape=1, size=2, width=0.1)+theme_bw()+xlab('FG poleward advance')+ ylab('Thermal midpoint (°C)')+
  theme(legend.position='none')

library(patchwork)
#png('C:/coral_fish/outputs/poleward_advance_plots.png',width =8, height =8 , units ="in", res =300)
p1+p2+p5+p3+p4+p6
dev.off()

###########################################################