### Loading packages ###
library(ggplot2)
library(FD)
library(ggrepel)
library(picante)
library(vegan)
library(geiger)
library(ape)
library(foreach)
library(apTreeshape)
library(MASS)
library(nlme)
library(AICcmodavg)
library(VGAM)
library(MuMIn)
library(scales)
library(car)
library(Hmisc)
###

### Load data ###
# Reading species richness for provinces with 3 and more species
abundance_t <- read.csv("Species_province.csv", header=T, row.names=1)

# Reading 14 traits data
traits <- read.csv("traits.csv", header=T, row.names=1)
str(traits) # check if binary traits are numeric as factors do not work for dbFD
traits[,2:14] <- sapply(traits[,2:14], as.numeric)

# Check that the order of species names are identical
names.traits <- rownames(traits)
names.abundance <- rownames(abundance_t)
identical(names.traits,names.abundance)
names.traits[!(names.traits %in% names.abundance)]
###


### calculate FD ###
abundance_FD <- t(abundance_t)

# Calculating FD
clams_FD <- dbFD(
  traits, abundance_FD, w.abun=TRUE, m="max" ,
  corr=c("cailliez"), clust.type="ward", 
  calc.FDiv = TRUE, calc.CWM = TRUE)

# Saving FD results
write.csv(clams_FD, file="output_clamFDresults.csv")
###


### PCoA plot ###
# Separate matrix for PCoA axes
pcoa_FD <- dbFD(
  traits, abundance_FD, w.abun=TRUE, m="max" ,
  corr=c("cailliez"), clust.type="ward", 
  calc.FDiv = TRUE, calc.CWM = TRUE, print.pco = TRUE)

# Plot PCoA
# Naming follows the loaded data
names_sp <- c("HP", "HH", "TG", "TMB", "TD", "TNO", "TC", "TS", "TM", "TR", "TSI", "TE")
ordiplot(pcoa_FD$x.axes)
text(pcoa_FD$x.axes, labels = names_sp, pos = 4)
###


### Alternative PCoA plot and vectors ###
x.dist <- gowdis(traits)
x.dist2 <- cailliez(x.dist)
x.pco <- dudi.pco(x.dist2, scannf=FALSE, full=TRUE)
li <- x.pco$li
points <- as.data.frame(li)

# percent difference explained by each Pco
sum(x.pco$eig[1:5]) / sum(x.pco$eig) #88.1%
sum(x.pco$eig[1:1]) / sum(x.pco$eig) #46.2%
sum(x.pco$eig[2:2]) / sum(x.pco$eig) #20.1%
sum(x.pco$eig[3:3]) / sum(x.pco$eig) #10.5%

colnames(traits) <- c("MD", "MSL", "FECD", "EGGSIZE", "FEMALEMAT", "GROWTH", "SYMMETRY", "HLSL", "UMBO", "RIB", "BYSGAPE", "SCUTE", "HYAL", "GUARDTEN")
ef <- envfit(x.pco, traits, permu = 999, na.rm=TRUE)
names <- row.names(traits)

plot(points[,1], points[,2], pch=21, xlab="PCoA 1 (46.2%)", ylab="PCoA 2 (20.1%)", bg=alpha("grey",0.5), xlim=c(-0.5,0.5), ylim=c(-0.6,0.5))
text(points, labels = names_sp, pos = 4)

plot(ef, p.max=1, col="blue", cex=.7) #spider
###


### Convex hull function for visualisation (optional) ###
Plot_ConvexHull<-function(xcoord, ycoord, lcolor){
  hpts <- chull(x = xcoord, y = ycoord)
  hpts <- c(hpts, hpts[1])
  polygon(xcoord[hpts], ycoord[hpts], col = lcolor)
} 

Plot_ConvexHull(xcoord = points[,1], ycoord = points[,2], lcolor = rgb(1,1,0.6,0.5)) #R1
Plot_ConvexHull(xcoord = points[c(4,5,7,9,10),1], ycoord = points[c(4,5,7,9,10),2], lcolor = rgb(1,0,0,0.5)) #R2
###


### Calculating null FD ###

clams_FD <- as.data.frame(clams_FD)
null.model.type <- "independentswap" #available methods: "independentswap", "frequency", "richness", "trialswap"

# create 999 permutations
nperm <- 999
FD_null <- matrix(NaN, nrow=nperm, ncol=length(clams_FD$FRic))
colnames(FD_null) <- rownames(clams_FD)
abund.null <- vector(mode="list", nperm)

# set seed for consistency
set.seed(20200525)

for(i in 1:nperm){ # loop for each permutation
  
  # Create permutated abundance matrix
  abund.i <- randomizeMatrix(abundance_FD, null.model = null.model.type)
  traits.i <- traits
  abund.null[[(i)]] <- abund.i
  
  # Record results
  res.i <- dbFD(
    traits.i, abund.i, w.abun=TRUE, 
    corr="cailliez", calc.FRic=TRUE, m="max",
    stand.FRic=FALSE, scale.RaoQ=FALSE, calc.FGR=FALSE,
    calc.CWM=FALSE, print.pco=FALSE
  )
  # Save result
  FD_null[i,] <- res.i$FRic
  
  rm(abund.i); rm(res.i)
  print(paste0(round(i/nperm*100), "% of permutations completed"))
}

FD_null_mean<-apply(FD_null, 2, mean)
FD_null_mean<-as.data.frame(FD_null_mean)
row.names(FD_null_mean)<-row.names(abundance_FD)
colnames(FD_null_mean)<- c("FD_null")

# wilcox test
clams_FRic <- matrix(NA, nrow=nrow(abundance_FD), ncol=1)
clams_FRic[,1] <- clams_FD$FRic
colnames(clams_FRic)<-c("FRic")
rownames(clams_FRic)<-rownames(abundance_FD)

wilcox_FD<-matrix(NA, nrow=nrow(abundance_FD), ncol=2)
colnames(wilcox_FD)<-c("p.value", "t.statistic")
rownames(wilcox_FD)<-rownames(abundance_FD)

for (i in 1:19){
  result <- wilcox.test(FD_null[,i], mu=clams_FRic[i,])
  wilcox_FD[i,1] <- result$p.value
  wilcox_FD[i,2] <- result$statistic
}
sum(wilcox_FD[,1]<0.05)
wilcox_FD[which(wilcox_FD[,1]>0.05),]

###

### Calculate PD ###

abundance_PD <- abundance_t
abunsp <- rownames(abundance_PD)
abunsp <- gsub(" ","_", abunsp)
rownames(abundance_PD) <- abunsp

# Species richness for provinces with 3 or more species
abundance_PD_all <- read.csv("Species_province.csv", header=T, row.names=1)
abunsp_all <- rownames(abundance_PD_all)
abunsp_all <- gsub(" ","_", abunsp_all)
rownames(abundance_PD_all) <- abunsp_all

# Reading nexus tree file
#use latest phylo tree
tree<-read.tree(file="Tridacninae.tre")

# Following labels would have to be changed since phylo tree is different
# Line below removes duplicate samples and reference genomes
tip<-c("TE-Ju40", "Trosewateri", "TL-33503", "WGT006", "WGT010", "WGT008", "WGT004", "WGT012", "TD", "TS", "FM", "AT", "CE")
newtree<-drop.tip(tree,tip)

newtree$tip.label[newtree$tip.label=="TE-Ju42"] <- "Tridacna_elongatissima"
newtree$tip.label[newtree$tip.label=="Tsi"] <- "Tridacna_squamosina"
newtree$tip.label[newtree$tip.label=="TL-33502"] <- "Tridacna_rosewateri"
newtree$tip.label[newtree$tip.label=="WGT005"] <- "Tridacna_maxima"
newtree$tip.label[newtree$tip.label=="WGT009"] <- "Tridacna_squamosa"
newtree$tip.label[newtree$tip.label=="WMT001"] <- "Tridacna_crocea"
newtree$tip.label[newtree$tip.label=="WGT007"] <- "Tridacna_noae"
newtree$tip.label[newtree$tip.label=="WGT003"] <- "Tridacna_derasa"
newtree$tip.label[newtree$tip.label=="WGT011"] <- "Tridacna_gigas"
newtree$tip.label[newtree$tip.label=="WGT013"] <- "Tridacna_mbalavuana"
newtree$tip.label[newtree$tip.label=="HP"] <- "Hippopus_porcellanus"
newtree$tip.label[newtree$tip.label=="HH"] <- "Hippopus_hippopus"

# Calculating phylogenetic diversity
clams_PD <- array(NA,dim=c(length(abundance_PD[1,]),1))
clams_PD <- as.data.frame(clams_PD)
rownames(clams_PD) <- colnames(abundance_PD)
colnames(clams_PD) <- c("PD")
clams_PD[,1] <-pd(t(abundance_PD_all),newtree)[,1]
shapiro.test(clams_PD[,1])
hist(clams_PD[,1])

write.csv(clams_PD, file="output_clamPDresults.csv")


# Calculating null PD 

# only generating one tree

null.model.type <- "independentswap" # available methods: "independentswap", "frequency", "richness", "trialswap"
nperm <- 999
PD_null <- matrix(NaN, nrow=nperm, ncol=nrow(clams_PD))
colnames(PD_null) <- colnames(abundance_PD)

# set seed for consistency
set.seed(20200526) 

for (i in seq(nperm)){
  abund.i <- randomizeMatrix(t(abundance_PD), null.model=null.model.type)
  tree <- newtree
  pd.i<- pd(abund.i, tree, include.root=FALSE)
  PD_null[i,] <- pd.i$PD
  rm(pd.i)
  print(paste0(round(i/nperm*100), "% of permutations completed"))
}

PD_null_mean <- apply(PD_null, 2, mean)
PD_null_mean <- as.matrix(PD_null_mean)
rownames(PD_null_mean) <- rownames(t(abundance_PD))
colnames(PD_null_mean) <- c("PD_null")

saveRDS(PD_null, file="data_PDnull.rds")

write.csv(PD_null, file= "clams_PDnull.csv")

# wilcox test
wilcox_PD<-matrix(NA, nrow=nrow(abundance_FD), ncol=2)
colnames(wilcox_PD)<-c("p.value", "t.statistic")
rownames(wilcox_PD)<-rownames(abundance_FD)

for (i in 1:19){
  result <- wilcox.test(PD_null[,i], mu=clams_PD[i,])
  wilcox_PD[i,1] <- result$p.value
  wilcox_PD[i,2] <- result$statistic
}
sum(wilcox_PD[,1]<0.05)

###

### Combining results ###

clams_results <- cbind(clams_FD[,1], clams_PD[,1], PD_null_mean[,1], clams_FD[,3], FD_null_mean[,1])
colnames(clams_results) <- c("SR", "PD", "null_PD", "FD", "null_FD")
rownames(clams_results) <- rownames(abundance_FD)
clams_results <- as.data.frame(clams_results)

#for null/observed heatmap
clams_results$PD_diff <- clams_results$PD-clams_results$null_PD
clams_results$FD_diff <- clams_results$FD-clams_results$null_FD

clams_results$PD_resid <- (clams_results$PD-clams_results$null_PD)/clams_results$PD
clams_results$FD_resid <- (clams_results$FD-clams_results$null_FD)/clams_results$FD

write.csv(clams_results, "total_clam_output.csv")
###

### Trait removal sensitivity analysis ###

# Create traits and weights for dataset with one trait removed
# Numbers represent the rows, columns in the .csv files
traits_nodepth    <- traits[,2:14]
traits_nolength   <- traits[,c(1,3:14)]
traits_nofecd     <- traits[,c(1:2,4:14)]
traits_noegg      <- traits[,c(1:3,5:14)]
traits_nofemal    <- traits[,c(1:4,6:14)]
traits_nogrowth   <- traits[,c(1:5,7:14)]
traits_nosymmetry <- traits[,c(1:6,8:14)]
traits_noHlSl     <- traits[,c(1:7,9:14)]
traits_noumbo     <- traits[,c(1:8,10:14)]
traits_norib      <- traits[,c(1:9,11:14)]
traits_nobysgape  <- traits[,c(1:10,12:14)]
traits_noscute    <- traits[,c(1:11,13:14)]
traits_nohyal     <- traits[,c(1:12,14:14)]
traits_noguardt   <- traits[,1:14]

# Calculate FD metrics, force m=2 to compare with original m="max"
res_nodepth    <- dbFD(traits_nodepth, abundance_FD, w.abun=TRUE, m=2, corr=c("cailliez"), clust.type="ward", calc.FDiv = TRUE)
res_nolength   <- dbFD(traits_nolength, abundance_FD, w.abun=TRUE, m=2, corr=c("cailliez"), clust.type="ward", calc.FDiv = TRUE)
res_nofecd     <- dbFD(traits_nofecd, abundance_FD, w.abun=TRUE, m=2, corr=c("cailliez"), clust.type="ward", calc.FDiv = TRUE)
res_noegg      <- dbFD(traits_noegg, abundance_FD, w.abun=TRUE, m=2, corr=c("cailliez"), clust.type="ward", calc.FDiv = TRUE)
res_nofemal    <- dbFD(traits_nofemal, abundance_FD, w.abun=TRUE, m=2, corr=c("cailliez"), clust.type="ward", calc.FDiv = TRUE)
res_nogrowth   <- dbFD(traits_nogrowth, abundance_FD, w.abun=TRUE, m=2, corr=c("cailliez"), clust.type="ward", calc.FDiv = TRUE)
res_nosymmetry <- dbFD(traits_nosymmetry, abundance_FD, w.abun=TRUE, m=2, corr=c("cailliez"), clust.type="ward", calc.FDiv = TRUE)
res_noHlSl     <- dbFD(traits_noHlSl, abundance_FD, w.abun=TRUE, m=2, corr=c("cailliez"), clust.type="ward", calc.FDiv = TRUE)
res_noumbo     <- dbFD(traits_noumbo, abundance_FD, w.abun=TRUE, m=2, corr=c("cailliez"), clust.type="ward", calc.FDiv = TRUE)
res_norib      <- dbFD(traits_norib, abundance_FD, w.abun=TRUE, m=2, corr=c("cailliez"), clust.type="ward", calc.FDiv = TRUE)
res_nobysgape  <- dbFD(traits_nobysgape, abundance_FD, w.abun=TRUE, m=2, corr=c("cailliez"), clust.type="ward", calc.FDiv = TRUE)
res_noscute    <- dbFD(traits_noscute, abundance_FD, w.abun=TRUE, m=2, corr=c("cailliez"), clust.type="ward", calc.FDiv = TRUE)
res_nohyal     <- dbFD(traits_nohyal, abundance_FD, w.abun=TRUE, m=2, corr=c("cailliez"), clust.type="ward", calc.FDiv = TRUE)
res_noguardt   <- dbFD(traits_noguardt, abundance_FD, w.abun=TRUE, m=2, corr=c("cailliez"), clust.type="ward", calc.FDiv = TRUE)

senslist <- list(as.data.frame(clams_FD), as.data.frame(res_nodepth), as.data.frame(res_nolength), as.data.frame(res_nofecd), as.data.frame(res_noegg), as.data.frame(res_nofemal), 
                 as.data.frame(res_nogrowth), as.data.frame(res_nosymmetry), as.data.frame(res_noHlSl), 
                 as.data.frame(res_noumbo), as.data.frame(res_norib), as.data.frame(res_nobysgape), 
                 as.data.frame(res_noscute), as.data.frame(res_nohyal), as.data.frame(res_noguardt))

sensFRic <- matrix(ncol=15, nrow=nrow(as.data.frame(clams_FD)))
sensnames <- c("observed","RM.DEPTH", "RM.LENGTH","RM.FECD", "RM.EGGSIZE", "RM.FEMALEMAT", "RM.GROWTH", "RM.SYMMETRY", "RM.HLSL", "RM.UMBO", "RM.RIB", "RM.BYSGAPE", "RM.SCUTE", "RM.HYAL", "RM.GUARDTEN")
colnames(sensFRic) <- sensnames

for (i in 1:length(senslist)){
  sensFRic[,i] <- senslist[[i]]$FRic
}

for (i in 1:(ncol(sensFRic)-1)) {
  print(cor.test(sensFRic[,1], sensFRic[,(i+1)], method="kendall"))
}

###

### Generalised linear model ###

###################################################################
# Linear Regression (FD~SR, PD~SR, PD~FD)##########################
###################################################################

shapiro.test(clams_results[,1]) #SR test for normality; p-value = 0.03157
shapiro.test(clams_results[,2]) #PD test for normality; p-value = 0.003564
shapiro.test(clams_results[,4]) #FD test for normality; p-value = 0.02067

clams_results <- as.data.frame(clams_results)
provinces <- rownames(clams_results)
clams_results <- cbind(clams_results, provinces)

realm <- read.csv("realm.csv", header=TRUE)
realm <- as.data.frame(realm)
rownames(realm) <- rownames(abundance_FD)
colnames(realm) <- c("provinces", "realm")
# realm <- subset(realm, rownames(realm)%in%row.names(clams_results)) 

clams_results <- merge(clams_results, realm, by="provinces", all.x=TRUE)
qqplot(clams_results$SR, clams_results$PD)
clams_results$lnPD <- sqrt(clams_results$PD)
qqplot(clams_results$SR, clams_results$lnPD)
shapiro.test(clams_results$PD)

# PD ~ SR
plot(clams_results$PD, clams_results$SR)
PD1 <- lm(PD~SR, data=clams_results)
PD2 <- gls(PD~SR, method="REML", data= clams_results)
vf1Fixed <- varFixed(~SR)
vf2 <- varPower(form=~SR)
PD2b <- gls(PD~SR, weights = vf1Fixed, method="REML", data= clams_results)
PD2c <- gls(PD~SR, weights = vf2, method="REML", data= clams_results)
PD3a <- lme(PD~SR, random = ~1|realm, method="REML", data= clams_results)
PD3b <- lme(PD~SR, random = ~1|provinces, method="REML", data= clams_results)
PD3c <- lme(PD~SR, random = ~1|realm/provinces, method="REML", data= clams_results)
AICc(PD1,PD3a,PD3b,PD3c)
summary(PD1) #0.7676

par(mfrow=c(2,2)) 
plot(PD1)
plot(resid(PD1)~SR,  data= clams_results)

#FRic ~ SR
plot(clams_results$FD, clams_results$SR)
FD1 <- lm(FD~SR, data=clams_results)
FD2 <- gls(FD~SR, method="REML", data= clams_results) # best model

FD3a <- lme(FD~SR, random = ~1|realm, method="REML", data= clams_results)
FD3b <- lme(FD~SR, random = ~1|provinces, method="REML", data= clams_results)
FD3c <- lme(FD~SR, random = ~1|realm/provinces, method="REML", data= clams_results)

AICc(FD1, FD3a, FD3b, FD3c)
par(mfrow=c(2,2))
plot(FD1)
summary(FD1) #0.7863

#FRic ~ PD
plot(clams_results$FD, clams_results$PD)
FP1 <- lm(FD~PD, data=clams_results)
FP2 <- gls(FD~PD, method="REML", data= clams_results) # best model

FP3a <- lme(FD~PD, random = ~1|realm, method="REML", data= clams_results)
FP3b <- lme(FD~PD, random = ~1|provinces, method="REML", data= clams_results)
FP3c <- lme(FD~PD, random = ~1|realm/provinces, method="REML", data= clams_results)

AICc(FP1, FP3a, FP3b, FP3c)
summary(FP1) #0.8179
par(mfrow=c(2,2))
plot(FP1)
###

### New outputs (start from here) ###
clams_results <- read.csv("total_clam_output.csv")

# glm_PD1 <- glm(PD~SR, data=clams_results, family = gaussian())
# summary(glm_PD1)
# plot(glm_PD1)

plot(clams_results$PD, clams_results$SR)

# Pearson's
clam_subset <- clams_results[,c("SR", "PD", "FD")]
cor(clam_subset[,c(1,2)], method="pearson") #SR~PD = 0.9298081
cor(clam_subset[,c(1,3)], method="pearson") #SR~FD = 0.9260012
cor(clam_subset[,c(2,3)], method="pearson") #PD~FD = 0.9727213

###

### Null plots ###

#PD transect level plot 
ggplot(data=clams_results, aes(SR, PD))+
  geom_point(shape=19, colour= "#999999")+labs(y=expression(PD["F"]))+
  geom_smooth(data=clams_results, aes(SR, null_PD), formula = y ~ x, method="loess", colour="black", se=F)+
  theme_bw()

#FRic transect level plot
ggplot(data=clams_results, aes(SR, FD))+
  geom_point(shape=19, colour= "#999999")+labs(y="FRic")+
  geom_smooth(data=clams_results, aes(SR, null_FD), formula = y ~ x, method='loess', colour="black", se=F)+
  theme_bw()

###


### Correlation plots ###

# Multiple plot function
#
# ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects)
# - cols:   Number of columns in layout
# - layout: A matrix specifying the layout. If present, 'cols' is ignored.
#
# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),
# then plot 1 will go in the upper left, 2 will go in the upper right, and
# 3 will go all the way across the bottom.
#
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)
  
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  
  numPlots = length(plots)
  
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  
  if (numPlots==1) {
    print(plots[[1]])
    
  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    
    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

p1 <- ggplot(clams_results, aes(SR, FD)) +
  geom_point(col="#999999") + theme_bw() + labs(y="FRic") +
  geom_smooth(formula = y ~ x, method="lm", col="blue", se=FALSE) +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14),
        axis.line = element_line(colour = "black"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.background = element_blank())

p2 <- ggplot(clams_results, aes(SR, PD)) +
  geom_point(col="#999999") + theme_bw()+labs(y=expression(PD["F"]))+
  geom_smooth(formula = y ~ x, method="lm", col="blue", se=FALSE)+
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14),
        axis.line = element_line(colour = "black"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.background = element_blank())

p3 <- ggplot(clams_results, aes(PD, FD)) +
  geom_point(col="#999999") + theme_bw()+labs(y=expression(FRic),x=expression(PD["F"]))+
  geom_smooth(formula = y ~ x, method="lm", col="blue", se=FALSE) +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14),
        axis.line = element_line(colour = "black"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.background = element_blank())

multiplot(p2,p1,p3, cols=3)
###

### Maps ###
#install.packages("PBSmapping")
library("PBSmapping")
library("ggplot2")
library("maps")
library("data.table")
library("stats")
library("ggsn")

# plot limits
xlim <- c(15, 240)
ylim <- c(-40, 40)
xlabel <- c("50", "100", "150", "-20")

worldmap <- map_data("world")
worldmap <- setNames(worldmap, c("X","Y","PID","POS","region","subregion"))
worldmap <- clipPolys(worldmap, xlim=xlim,ylim=ylim, keepExtra=TRUE)

# All provinces with nbsp
provinces_SR <- read.csv ("provinces_SR.csv")
# provinces.stats from Meow_map.R

#test reverse contrast

PD_resid_province_ggplot_col <- ggplot(provinces.stats) + 
  theme(panel.background = element_rect(fill="grey85"),
        panel.grid.major = element_line(colour="grey90"),
        legend.position="bottom") +
  aes(long.recenter,lat,group=group) + 
  geom_polygon(aes(fill=PD_resid)) +
  scale_fill_gradientn(colors = c("#3182bd", "white", "red"), values = scales::rescale(c(-0.518626614, 0, 0.231064967)), guide="colourbar") +
  coord_equal() 

FD_resid_province_ggplot_col <- ggplot(provinces.stats) + 
  theme(panel.background = element_rect(fill="grey85"),
        panel.grid.major = element_line(colour="grey90"),
        legend.position="bottom") +
  aes(long.recenter,lat,group=group) + 
  geom_polygon(aes(fill=FD_resid)) +
  scale_fill_gradientn(colors = c("#3182bd", "white", "red"), values = scales::rescale(c(-5.002282219, 0, 0.473492253)), guide="colourbar") +
  coord_equal()

# Clam FD residuals difference by Province
p13 <- FD_resid_province_ggplot_col +
  coord_map(xlim=xlim,ylim=ylim) +
  geom_polygon(data=worldmap,aes(X,Y,group=PID),
               fill = "grey95",color="black", size=0.2) +
  labs(y="",x="") +
  scale_x_continuous(breaks = c(50, 100, 150, 200), labels = xlabel)

# Clam PD residuals difference by Province
p14 <- PD_resid_province_ggplot_col +
  coord_map(xlim=xlim,ylim=ylim) +
  geom_polygon(data=worldmap,aes(X,Y,group=PID),
               fill = "grey95",color="black", size=0.2) +
  labs(y="",x="") +
  scale_x_continuous(breaks = c(50, 100, 150, 200), labels = xlabel)

# adding scale bars 
anchor <- c(x=231, y=-35)
pd <- p13 + scalebar(provinces.stats, dist = 2000, dist_unit = "km", transform = TRUE, model = "WGS84", location ="bottomright", anchor = anchor, st.size = 4, st.dist = 0.03, border.size = 0.5)
pe <- p14 + scalebar(provinces.stats, dist = 2000, dist_unit = "km", transform = TRUE, model = "WGS84", location ="bottomright", anchor = anchor, st.size = 4, st.dist = 0.03, border.size = 0.5)
#pf <- p15 + scalebar(provinces.stats, dist = 2000, dist_unit = "km", transform = TRUE, model = "WGS84", location ="bottomright", anchor = anchor, st.size = 4, st.dist = 0.03, border.size = 0.5)

# PD and FD resid plots
multiplot(pe, pd)
###