
library(ggplot2)
library(ggridges)
library(grid)
library(data.table)
library(reshape2)

library(Hmsc)
library(reshape2)
library(abind)
library(tidyr)
library(dplyr)
#library(berryFunctions)
#library(lemon)
#library(cowplot)
library(patchwork)
library(ggh4x)
library(ggpubr)

library(vioplot)

setwd("D:/Helsinki/Moth_project/ResultsAug23/RDATA/3.ModelRun")
load('VPbwFIN_PA01408.RDATA')
load('VPcondAFIN100cdd_PA_1408.RDATA')
dim(VP.bw$vals)

## a
VP.a.M1H = VP.a$vals#[,as.numeric(rownames(M1H)),]
dim(VP.a.M1H)
M1Ha.vp = apply(VP.a.M1H, c(1,3), mean)
dim(M1Ha.vp)

M1Ha.low = apply(VP.a.M1H, c(1,3), function(x) quantile(x, prob=0.025))
M1Ha.up = apply(VP.a.M1H, c(1,3), function(x) quantile(x, prob=0.975))



## bw
VP.bw.M1H = VP.bw$vals#[,as.numeric(rownames(M1H)),,]
dim(VP.bw.M1H)
M1Hbw.vp = apply(VP.bw.M1H, c(1,3,4), mean)
dim(M1Hbw.vp)

M1Hbw.low = apply(VP.bw.M1H, c(1,3,4), function(x) quantile(x, prob=0.025))
M1Hbw.up = apply(VP.bw.M1H, c(1,3,4), function(x) quantile(x, prob=0.975))

# mean over species and MCMC samples
Meanbw.vp1 = apply(VP.bw.M1H, c(1,3), mean)


## cond
VP.cond.M1H = VP.cond$vals#[,as.numeric(rownames(M1H)),,]
dim(VP.cond.M1H)
M1Hcond.vp = apply(VP.cond.M1H, c(1,3,4), mean)
dim(M1Hcond.vp)

M1Hcond.low = apply(VP.cond.M1H, c(1,3,4), function(x) quantile(x, prob=0.025))
M1Hcond.up = apply(VP.cond.M1H, c(1,3,4), function(x) quantile(x, prob=0.975))

# mean over species and MCMC samples
Meancond.vp1 = apply(VP.cond.M1H, c(1,3), mean)

{## transform data
  
  ### cond
  #### mean
  dim(M1Hcond.vp)
  rownames(M1Hcond.vp)
  
  M1H.valcond = M1Hcond.vp
  M1H.valcond <- aperm(M1Hcond.vp, c(2,1,3))
  rownames(M1H.valcond) = VP.cond$G
  
  valcondM1H = as.data.frame(apply(M1H.valcond, 2, c))
  
  valcondM1H$Cond = as.vector(rep(rownames(M1H.valcond), dim(M1Hcond.vp)[3]))
  
  valcond.M1H_long <- reshape2::melt(valcondM1H)
  colnames(valcond.M1H_long) = c('Cond', 'Var', 'VP')
  valcond.M1H_long$group = rep('Group', dim(valcond.M1H_long)[1])
  valcond.M1H_long$measure = rep('mean', dim(valcond.M1H_long)[1])
  
  #### up
  dim(M1Hcond.up)
  rownames(M1Hcond.up)
  
  M1H.valcond = M1Hcond.up
  M1H.valcondup <- aperm(M1Hcond.up, c(2,1,3))
  rownames(M1H.valcondup) = VP.cond$G
  
  valcondM1H.up = as.data.frame(apply(M1H.valcondup, 2, c))
  
  valcondM1H.up$Cond = as.vector(rep(rownames(M1H.valcondup), dim(M1Hcond.up)[3]))
  
  valcond.M1H.up_long <- reshape2::melt(valcondM1H.up)
  colnames(valcond.M1H.up_long) = c('Cond', 'Var', 'VP')
  valcond.M1H.up_long$group = rep('Group', dim(valcond.M1H.up_long)[1])
  valcond.M1H.up_long$measure = rep('CI_up', dim(valcond.M1H.up_long)[1])
  
  #### low
  dim(M1Hcond.low)
  rownames(M1Hcond.low)
  
  M1H.valcond = M1Hcond.low
  M1H.valcondlow <- aperm(M1Hcond.low, c(2,1,3))
  rownames(M1H.valcondlow) = VP.cond$G
  
  valcondM1H.low = as.data.frame(apply(M1H.valcondlow, 2, c))
  
  valcondM1H.low$Cond = as.vector(rep(rownames(M1H.valcondlow), dim(M1Hcond.low)[3]))
  
  valcond.M1H.low_long <- reshape2::melt(valcondM1H.low)
  colnames(valcond.M1H.low_long) = c('Cond', 'Var', 'VP')
  valcond.M1H.low_long$group = rep('Group', dim(valcond.M1H.low_long)[1])
  valcond.M1H.low_long$measure = rep('CI_low', dim(valcond.M1H.low_long)[1])
  
  
  ### bw
  #### mean
  dim(M1Ha.vp)
  
  M1H.valbw.tmp <- aperm(M1Hbw.vp, c(1,3,2))
  dim(M1H.valbw.tmp)
  M1Hbw.vp2 = abind(M1H.valbw.tmp, along = 3)
  
  M1H.valbw <- aperm(M1Hbw.vp2, c(3,1,2))
  rownames(M1H.valbw) = c('BTW', 'WITHIN')
  
  valbwM1H = as.data.frame(apply(M1H.valbw, 2, c))
  
  valbwM1H$Cond = as.vector(rep(rownames(M1H.valbw), dim(M1Hbw.vp)[3]))
  
  valbw.M1H_long <- reshape2::melt(valbwM1H)
  colnames(valbw.M1H_long) = c('Cond', 'Var', 'VP')
  valbw.M1H_long$group = rep('BW', dim(valbw.M1H_long)[1])
  valbw.M1H_long$measure = rep('mean', dim(valbw.M1H_long)[1])
  
  #### up
  dim(M1Ha.up)
  
  M1H.valbwup.tmp <- aperm(M1Hbw.up, c(1,3,2))
  dim(M1H.valbwup.tmp)
  M1Hbw.up2 = abind(M1H.valbwup.tmp,along = 3)
  
  M1H.valbwup <- aperm(M1Hbw.up2, c(3,1,2))
  rownames(M1H.valbwup) = c('BTW', 'WITHIN')
  
  valbwM1H.up = as.data.frame(apply(M1H.valbwup, 2, c))
  
  valbwM1H.up$Cond = as.vector(rep(rownames(M1H.valbwup), dim(M1Hbw.up)[3]))
  
  valbw.M1H.up_long <- reshape2::melt(valbwM1H.up)
  colnames(valbw.M1H.up_long) = c('Cond', 'Var', 'VP')
  valbw.M1H.up_long$group = rep('BW', dim(valbw.M1H.up_long)[1])
  valbw.M1H.up_long$measure = rep('CI_up', dim(valbw.M1H.up_long)[1])
  
  #### low
  dim(M1Ha.low)
  
  M1H.valbwlow.tmp <- aperm(M1Hbw.low, c(1,3,2))
  dim(M1H.valbwlow.tmp)
  M1Hbw.low2 = abind(M1H.valbwlow.tmp, along = 3)
  dim(M1Hbw.low2)
  
  M1H.valbwlow <- aperm(M1Hbw.low2, c(3,1,2))
  rownames(M1H.valbwlow) = c('BTW', 'WITHIN')
  
  valbwM1H.low = as.data.frame(apply(M1H.valbwlow, 2, c))
  
  valbwM1H.low$Cond = as.vector(rep(rownames(M1H.valbwlow), dim(M1Hbw.low)[3]))
  
  valbw.M1H.low_long <- reshape2::melt(valbwM1H.low)
  colnames(valbw.M1H.low_long) = c('Cond', 'Var', 'VP')
  valbw.M1H.low_long$group = rep('BW', dim(valbw.M1H.low_long)[1])
  valbw.M1H.low_long$measure = rep('CI_low', dim(valbw.M1H.low_long)[1])
  
  
  # combine
  VP.all1 = rbind(valcond.M1H_long, valbw.M1H_long, 
                  valcond.M1H.up_long, valbw.M1H.up_long,
                  valcond.M1H.low_long,valbw.M1H.low_long)
  
  VP.all1 = rbind(valcond.M1H_long, valbw.M1H_long)
  head(VP.all1)
  
  VP.all1$Cond = factor(VP.all1$Cond, levels=c('Heterogeneous', 'Forest', 'Wetland - vegetation', 'Water', 'BTW', 'WITHIN'))
}

VPPA.alB_reg = VP.all1

###  AB
load('VPbwFIN_AB01608.RDATA')
load('VPcondAFIN100cdd_AB_1408.RDATA')
dim(VP.bw$vals)

## a
VP.a.M1H = VP.a$vals#[,as.numeric(rownames(M1H)),]
dim(VP.a.M1H)
M1Ha.vp = apply(VP.a.M1H, c(1,3), mean)
dim(M1Ha.vp)

M1Ha.low = apply(VP.a.M1H, c(1,3), function(x) quantile(x, prob=0.025))
M1Ha.up = apply(VP.a.M1H, c(1,3), function(x) quantile(x, prob=0.975))



## bw
VP.bw.M1H = VP.bw$vals#[,as.numeric(rownames(M1H)),,]
dim(VP.bw.M1H)
M1Hbw.vp = apply(VP.bw.M1H, c(1,3,4), mean)
dim(M1Hbw.vp)

M1Hbw.low = apply(VP.bw.M1H, c(1,3,4), function(x) quantile(x, prob=0.025))
M1Hbw.up = apply(VP.bw.M1H, c(1,3,4), function(x) quantile(x, prob=0.975))

# mean over species and MCMC samples
Meanbw.vp1 = apply(VP.bw.M1H, c(1,3), mean)


## cond
VP.cond.M1H = VP.cond$vals#[,as.numeric(rownames(M1H)),,]
dim(VP.cond.M1H)
M1Hcond.vp = apply(VP.cond.M1H, c(1,3,4), mean)
dim(M1Hcond.vp)

M1Hcond.low = apply(VP.cond.M1H, c(1,3,4), function(x) quantile(x, prob=0.025))
M1Hcond.up = apply(VP.cond.M1H, c(1,3,4), function(x) quantile(x, prob=0.975))

# mean over species and MCMC samples
Meancond.vp1 = apply(VP.cond.M1H, c(1,3), mean)

{## transform data
  
  ### cond
  #### mean
  dim(M1Hcond.vp)
  rownames(M1Hcond.vp)
  
  M1H.valcond = M1Hcond.vp
  M1H.valcond <- aperm(M1Hcond.vp, c(2,1,3))
  rownames(M1H.valcond) = VP.cond$G
  
  valcondM1H = as.data.frame(apply(M1H.valcond, 2, c))
  
  valcondM1H$Cond = as.vector(rep(rownames(M1H.valcond), dim(M1Hcond.vp)[3]))
  
  valcond.M1H_long <- reshape2::melt(valcondM1H)
  colnames(valcond.M1H_long) = c('Cond', 'Var', 'VP')
  valcond.M1H_long$group = rep('Group', dim(valcond.M1H_long)[1])
  valcond.M1H_long$measure = rep('mean', dim(valcond.M1H_long)[1])
  
  #### up
  dim(M1Hcond.up)
  rownames(M1Hcond.up)
  
  M1H.valcond = M1Hcond.up
  M1H.valcondup <- aperm(M1Hcond.up, c(2,1,3))
  rownames(M1H.valcondup) = VP.cond$G
  
  valcondM1H.up = as.data.frame(apply(M1H.valcondup, 2, c))
  
  valcondM1H.up$Cond = as.vector(rep(rownames(M1H.valcondup), dim(M1Hcond.up)[3]))
  
  valcond.M1H.up_long <- reshape2::melt(valcondM1H.up)
  colnames(valcond.M1H.up_long) = c('Cond', 'Var', 'VP')
  valcond.M1H.up_long$group = rep('Group', dim(valcond.M1H.up_long)[1])
  valcond.M1H.up_long$measure = rep('CI_up', dim(valcond.M1H.up_long)[1])
  
  #### low
  dim(M1Hcond.low)
  rownames(M1Hcond.low)
  
  M1H.valcond = M1Hcond.low
  M1H.valcondlow <- aperm(M1Hcond.low, c(2,1,3))
  rownames(M1H.valcondlow) = VP.cond$G
  
  valcondM1H.low = as.data.frame(apply(M1H.valcondlow, 2, c))
  
  valcondM1H.low$Cond = as.vector(rep(rownames(M1H.valcondlow), dim(M1Hcond.low)[3]))
  
  valcond.M1H.low_long <- reshape2::melt(valcondM1H.low)
  colnames(valcond.M1H.low_long) = c('Cond', 'Var', 'VP')
  valcond.M1H.low_long$group = rep('Group', dim(valcond.M1H.low_long)[1])
  valcond.M1H.low_long$measure = rep('CI_low', dim(valcond.M1H.low_long)[1])
  
  
  ### bw
  #### mean
  dim(M1Ha.vp)
  
  M1H.valbw.tmp <- aperm(M1Hbw.vp, c(1,3,2))
  dim(M1H.valbw.tmp)
  M1Hbw.vp2 = abind(M1H.valbw.tmp,along = 3)
  
  M1H.valbw <- aperm(M1Hbw.vp2, c(3,1,2))
  rownames(M1H.valbw) = c('BTW', 'WITHIN')
  
  valbwM1H = as.data.frame(apply(M1H.valbw, 2, c))
  
  valbwM1H$Cond = as.vector(rep(rownames(M1H.valbw), dim(M1Hbw.vp)[3]))
  
  valbw.M1H_long <- reshape2::melt(valbwM1H)
  colnames(valbw.M1H_long) = c('Cond', 'Var', 'VP')
  valbw.M1H_long$group = rep('BW', dim(valbw.M1H_long)[1])
  valbw.M1H_long$measure = rep('mean', dim(valbw.M1H_long)[1])
  
  #### up
  dim(M1Ha.up)
  
  M1H.valbwup.tmp <- aperm(M1Hbw.up, c(1,3,2))
  dim(M1H.valbwup.tmp)
  M1Hbw.up2 = abind(M1H.valbwup.tmp, along = 3)
  
  M1H.valbwup <- aperm(M1Hbw.up2, c(3,1,2))
  rownames(M1H.valbwup) = c('BTW', 'WITHIN')
  
  valbwM1H.up = as.data.frame(apply(M1H.valbwup, 2, c))
  
  valbwM1H.up$Cond = as.vector(rep(rownames(M1H.valbwup), dim(M1Hbw.up)[3]))
  
  valbw.M1H.up_long <- reshape2::melt(valbwM1H.up)
  colnames(valbw.M1H.up_long) = c('Cond', 'Var', 'VP')
  valbw.M1H.up_long$group = rep('BW', dim(valbw.M1H.up_long)[1])
  valbw.M1H.up_long$measure = rep('CI_up', dim(valbw.M1H.up_long)[1])
  
  #### low
  dim(M1Ha.low)
  
  M1H.valbwlow.tmp <- aperm(M1Hbw.low, c(1,3,2))
  dim(M1H.valbwlow.tmp)
  M1Hbw.low2 = abind(M1H.valbwlow.tmp, along = 3)
  dim(M1Hbw.low2)
  
  M1H.valbwlow <- aperm(M1Hbw.low2, c(3,1,2))
  rownames(M1H.valbwlow) = c('BTW', 'WITHIN')
  
  valbwM1H.low = as.data.frame(apply(M1H.valbwlow, 2, c))
  
  valbwM1H.low$Cond = as.vector(rep(rownames(M1H.valbwlow), dim(M1Hbw.low)[3]))
  
  valbw.M1H.low_long <- reshape2::melt(valbwM1H.low)
  colnames(valbw.M1H.low_long) = c('Cond', 'Var', 'VP')
  valbw.M1H.low_long$group = rep('BW', dim(valbw.M1H.low_long)[1])
  valbw.M1H.low_long$measure = rep('CI_low', dim(valbw.M1H.low_long)[1])
  
  
  # combine
  VP.all1 = rbind(valcond.M1H_long, valbw.M1H_long, 
                  valcond.M1H.up_long, valbw.M1H.up_long,
                  valcond.M1H.low_long,valbw.M1H.low_long)
  
  VP.all1 = rbind(valcond.M1H_long, valbw.M1H_long)
  head(VP.all1)
  
  VP.all1$Cond = factor(VP.all1$Cond, levels=c('Heterogeneous', 'Forest', 'Wetland - vegetation', 'Water', 'BTW', 'WITHIN'))
}

VPAB.alB_reg = VP.all1




############################################## VP from model
setwd("D:/Helsinki/Moth_project/ResultsAug23/RDATA/1.EnvironmentalInfo")
load( 'clusterdata2.RDATA')

head(mapData3)
mapData3 = mapData3[order(mapData3$X.x),]
# Need for all level of the sampling design!!!
DF = mapData3[, c(1, 20)]
colnames(DF)[1] = 'sample'
DF[,1] = gsub(" ", " _ ", DF[,1])

head(DF)

colors =c("#EE6677", "#228833", 
          "#66CCEE", "#CCBB44",
          '#BBBBBB', 'black')
colBR = data.frame(Cond=c('Heterogeneous', 'Forest', 'Wetland - vegetation', 'Water', 'BTW', 'WITHIN'), col=colors)



VPPA.alB_reg.mean = VPPA.alB_reg[VPPA.alB_reg$measure == 'mean',]
VPAB.alB_reg.mean = VPAB.alB_reg[VPAB.alB_reg$measure == 'mean',]


# year
VPreg_y = cbind(VPPA.alB_reg.mean, VPAB.alB_reg.mean$VP)
colnames(VPreg_y)[6] = 'VPAB'
head(VPreg_y)

colnames(colBR)[1] = 'Cond'
VPreg_y.col = merge(VPreg_y, colBR, by='Cond')

           
################################################# Overall



library(maps)
library(mapdata)

library(sp)
#library(rgdal)
library(raster)
library(viridis)

Finmap = map('worldHires','Finland')
Fin2 <- maps::map('worldHires','Finland', plot=FALSE, fill=TRUE)

zones=matrix(c(0,3,11,6,12, 1,2,4,5,7, 8,9,10,13,14), ncol=5, byrow=TRUE)
layout(zones, widths=c(3,5,1,5,1), heights = c(1,4,1))

# 
par(mai=c(0.05, 0.05, 0.051, 0.05))
maps::map(Fin2, fill=TRUE, col='gray60')

mapData3$Cond = mapData3$col
FIN.spdf <- SpatialPointsDataFrame(mapData3[,c(2,3)], mapData3)
proj4string(FIN.spdf) <-CRS("+proj=utm +zone=35 +ellps=GRS80 +units=m +no_defs")

newFin<-spTransform(FIN.spdf, CRS("+init=epsg:4326"))

siteCol = merge(mapData3, colBR, by='Cond')
siteCol = siteCol[match(Site, siteCol$Site),]

points(newFin, col=siteCol$col.y, pch=16, cex=1.6)

# ------------------------------------------------------------------- Old plotting
#par(mai=c(0.78, 0.78, 1, 0.45), xpd=F)


# ## VPyear PA
# 
# par(mar=c(5,5,1,1))
# #VPreg_y.Clim = VPreg_y.col[VPreg_y.col$Var == 'Climate',]
# plot(VPreg_y.col[VPreg_y.col$Var=='Climate',]$VP ~ VPreg_y.col[VPreg_y.col$Var=='Habitat',]$VP, 
#      col=VPreg_y.col[VPreg_y.col$Var=='Climate',]$col,# main = 'PA Year',
#      xlab='VP-Habitat', ylab='VP-Climate', cex.axis = 1.5, cex.lab = 2, cex.main=2,
#      xlim = c(0, 0.55), ylim=c(0, 0.55))
# abline(coef = c(0,1), col='black')
# 
# par(mar=c(0,5,1,1))
# 
# VPreg_y.H = VPreg_y.col[VPreg_y.col$Var == 'Habitat',]
# VPreg_y.H.BW = VPreg_y.H[VPreg_y.H$Cond %in% c('BTW', 'WITHIN'),]
# vioplot(VPreg_y.H[VPreg_y.H$Cond == 'BTW',]$VP, bty = 'n', axes=F,
#         ylim = c(0, 0.5), xlim=c(1,1.5), space=0, horiz=TRUE, col=alpha('yellow', 0.6),
#         horizontal=T, plotCentre = "line", side = 'right',
#         yaxt='n', xaxt='n', ann=FALSE)
# vioplot(VPreg_y.H[VPreg_y.H$Cond == 'WITHIN',]$VP, bty = 'n', axes=F,
#         ylim = c(0, 0.5), space=0, horiz=TRUE, col=alpha('darkorange', 0.6),
#         horizontal=T, plotCentre = "line", side = 'right', add=T)
# 
# 
# par(mar=c(5,0,1,1))
# VPreg_y.C = VPreg_y.col[VPreg_y.col$Var == 'Climate',]
# VPreg_y.C.BW = VPreg_y.C[VPreg_y.C$Cond %in% c('BTW', 'WITHIN'),]
# vioplot(VPreg_y.C.BW[VPreg_y.C.BW$Cond == 'BTW',]$VP, bty = 'n', axes=F,
#         ylim = c(0, 0.5), xlim=c(1,1.5), space=0, horiz=TRUE, col=alpha('yellow', 0.6),
#         plotCentre = "line", side = 'right', xaxt='n', ann=FALSE, yaxt='n')
# vioplot(VPreg_y.C.BW[VPreg_y.C.BW$Cond == 'WITHIN',]$VP, bty = 'n', axes=F,
#         ylim = c(0, 0.5), space=0, horiz=TRUE, col=alpha('darkorange', 0.6),
#         plotCentre = "line", side = 'right', add=T)
# 
# 
# 
# ## VPyear AB
# 
# par(mar=c(5,5,1,1))
# #VPreg_y.Clim = VPreg_y.col[VPreg_y.col$Var == 'Climate',]
# plot(VPreg_y.col[VPreg_y.col$Var=='Climate',]$VPAB ~ VPreg_y.col[VPreg_y.col$Var=='Habitat',]$VPAB, 
#      col=VPreg_y.col[VPreg_y.col$Var=='Climate',]$col, #main = 'AB Year',
#      xlab='VP-Habitat', ylab='VP-Climate', cex.axis = 1.5, cex.lab = 2, cex.main=2,
#      xlim = c(0, 0.5), ylim=c(0, 0.5))
# abline(coef = c(0,1), col='black')
# 
# par(mar=c(0,5,1,1))
# 
# VPreg_y.H = VPreg_y.col[VPreg_y.col$Var == 'Habitat',]
# VPreg_y.H.BW = VPreg_y.H[VPreg_y.H$Cond %in% c('BTW', 'WITHIN'),]
# vioplot(VPreg_y.H[VPreg_y.H$Cond == 'BTW',]$VPAB, bty = 'n', axes=F,
#         ylim = c(0, 0.5), xlim=c(1,1.5), space=0, horiz=TRUE, col=alpha('yellow', 0.6),
#         horizontal=T, plotCentre = "line", side = 'right',
#         yaxt='n', xaxt='n', ann=FALSE)
# vioplot(VPreg_y.H[VPreg_y.H$Cond == 'WITHIN',]$VPAB, bty = 'n', axes=F,
#         ylim = c(0, 0.5), space=0, horiz=TRUE, col=alpha('darkorange', 0.6),
#         horizontal=T, plotCentre = "line", side = 'right', add=T)
# 
# 
# par(mar=c(5,0,1,1))
# VPreg_y.C = VPreg_y.col[VPreg_y.col$Var == 'Climate',]
# VPreg_y.C.BW = VPreg_y.C[VPreg_y.C$Cond %in% c('BTW', 'WITHIN'),]
# vioplot(VPreg_y.C.BW[VPreg_y.C.BW$Cond == 'BTW',]$VPAB, bty = 'n', axes=F,
#         ylim = c(0, 0.5), xlim=c(1,1.5), space=0, horiz=TRUE, col=alpha('yellow', 0.6),
#         plotCentre = "line", side = 'right', xaxt='n', ann=FALSE, yaxt='n')
# vioplot(VPreg_y.C.BW[VPreg_y.C.BW$Cond == 'WITHIN',]$VPAB, bty = 'n', axes=F,
#         ylim = c(0, 0.5), space=0, horiz=TRUE, col=alpha('darkorange', 0.6),
#         plotCentre = "line", side = 'right', add=T)
# 
# 
# 
# # legend upper square
# par(mar=c(1,0,1,1))
# plot(0,0, bty='n', col='white',axes=F, yaxt='n', xaxt='n', ylab='', xlab='')
# legend(legend=c('BTW', 'WITHIN'), 
#        'top', pch=22, pt.cex=5,
#        pt.bg=colors[5:6], bty = "n", cex = 1.3,
#        col=colors[5:6])
# 
# # legend lower square
# plot(0,0, bty='n', col='white',axes=F, yaxt='n', xaxt='n', ylab='', xlab='')
# legend(legend=c('Heterog', 'Forest', 'WaterVG', 'Urban'), 
#        'center', pch=22, pt.cex=5, ncol=2,
#        pt.bg=colors[1:4], bty = "n", cex = 1.3,
#        col=colors[1:4])
# 
# 
# 
# # Plots names
# plot(0,0, bty='n', col='white',axes=F, yaxt='n', xaxt='n', ylab='', xlab='')
# title("Year model - AB", line = -6, outer = T, cex.main = 2, adj=0.78)
# plot(0,0, bty='n', col='white',axes=F, yaxt='n', xaxt='n', ylab='', xlab='')
# title("Year model - PA", line = -6, outer = T, cex.main = 2, adj=0.25)
# 
# 
# plot(0,0, bty='n', col='white',axes=F, yaxt='n', xaxt='n', ylab='', xlab='')
# 


#.................................................................................................






VPPA.alB_reg.mean = VPPA.alB_reg#[VPPA.alB_reg$Cond %in% c('Heterogeneous', 'Forest', 
                                                       # 'Water', 'Wetland - vegetation'),]
VPAB.alB_reg.mean = VPAB.alB_reg#[VPAB.alB_reg$Cond %in% c('Heterogeneous', 'Forest', 
                                                       #   'Water', 'Wetland - vegetation'),]



####
library(ggplot2)
library(ggExtra)
library(cowplot)

# With ggplot

head(VPPA.alB_reg.mean)

VPPA.alB_reg.mean = VPPA.alB_reg.mean[VPPA.alB_reg.mean$measure %in% 'mean',]

VPreg_yPA.hab = VPPA.alB_reg.mean[VPPA.alB_reg.mean$Var %in% c('Habitat'),]
VPreg_yPA.clim = VPPA.alB_reg.mean[VPPA.alB_reg.mean$Var %in% c('Climate'),]

VPreg_yPA = cbind(VPreg_yPA.clim, VPreg_yPA.hab$VP)
colnames(VPreg_yPA)[c(3, 6)] = c('Climate', 'Habitat')


gPA = ggplot(VPreg_yPA, aes(x = Habitat , y = Climate, fill = Cond)) +
  #geom_point()+
  stat_density_2d(aes(alpha = ..level.., fill = Cond), geom = "polygon", 
                  position = "identity", h=0.05, n=50, bins=50, show.legend = FALSE) +
  #scale_color_manual(values = c("forestgreen", "brown2", "azure3", "darkcyan"))+
  scale_fill_manual(values = c('Heterogeneous'="#EE6677", 'Forest'="#228833", 
                               'Water'="#66CCEE", 'Wetland - vegetation'="#CCBB44",
                               'BTW'='#BBBBBB', 'WITHIN'='black'))+
  theme_pubr(legend = c(0.25,0.9)) + 
  ylim(0.05, 0.78) + xlim(0.05,0.78) +
  labs(
    x = "VP habitat",
    y = "VP climate",
    fill = "Cluster"
  ) + theme(plot.margin = unit(c(0, 0, 0, 0), "cm"),
            legend.title = element_text(size=14), 
            legend.text = element_text(size=14),
            axis.title.x = element_text(size = 16),
            axis.text.x = element_text(size = 14),
            axis.text.y = element_text(size = 14),
            axis.title.y = element_text(size = 16)) +
  geom_abline(intercept = 0, slope = 1) + guides(alpha='none') 


VPreg_yPA.col.sub = VPPA.alB_reg.mean[VPPA.alB_reg.mean$Cond %in% c('BTW', 'WITHIN'),]
VPreg_yPA.col.hab = VPreg_yPA.col.sub[VPreg_yPA.col.sub$Var %in% c('Habitat'),]
VPreg_yPA.col.clim = VPreg_yPA.col.sub[VPreg_yPA.col.sub$Var %in% c('Climate'),]

xmarg <- ggplot(VPreg_yPA.col.hab, aes(x = VP, y=..density.., fill = Cond)) + 
  geom_density(alpha = 0.4)  + xlim(0.05, 0.78) +
  scale_fill_manual(values = c('#BBBBBB','black')) + 
  theme_void() + theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) + 
  theme(legend.position = "none")+ ggtitle('PA model') +
  theme(plot.title = element_text(hjust = 1.2, vjust =  -15, size=14)) 


ymarg <- ggplot(VPreg_yPA.col.clim) + 
  geom_density(alpha = 0.4, aes(x = VP, fill = Cond)) + xlim(0.05, 0.78) +
  scale_fill_manual(values = c('#BBBBBB','black')) + 
  theme_void()  + theme(plot.margin = unit(c(0, 0, 0, 0), "cm"))+ 
  theme(legend.position = "none") + 
  coord_flip()





VPreg_yAB.hab = VPAB.alB_reg.mean[VPAB.alB_reg.mean$Var %in% c('Habitat'),]
VPreg_yAB.clim = VPAB.alB_reg.mean[VPAB.alB_reg.mean$Var %in% c('Climate'),]

VPreg_yAB = cbind(VPreg_yAB.clim, VPreg_yAB.hab$VP)
colnames(VPreg_yAB)[c(3, 6)] = c('Climate', 'Habitat')

gAB = ggplot(VPreg_yAB, aes(x = Habitat , y = Climate, fill = Cond)) +
  #geom_point()+
  stat_density_2d(aes(alpha = ..level.., fill = Cond), geom = "polygon", 
                  position = "identity", h=0.05, n=50, bins=50) +
  #scale_color_manual(values = c("forestgreen", "brown2", "azure3", "darkcyan"))+
  scale_fill_manual(values = c('Heterogeneous'="#EE6677", 'Forest'="#228833", 
                               'Water'="#66CCEE", 'Wetland - vegetation'="#CCBB44",
                               'BTW'='#BBBBBB', 'WITHIN'='black'))+
  theme_pubr(legend = c(0.85,0.9)) + 
  ylim(0.05, 0.78) + xlim(0.05,0.78) +
  labs(
    x = "VP habitat",
    y = "VP climate",
    fill = "Cluster"
  ) + theme(plot.margin = unit(c(0, 0, 0, 0), "cm"),
            legend.title = element_text(size=14), 
            legend.text = element_text(size=14),
            axis.title.x = element_text(size = 16),
            axis.text.x = element_text(size = 14),
            axis.text.y = element_text(size = 14),
            axis.title.y = element_text(size = 16)) +
  geom_abline(intercept = 0, slope = 1) + guides(alpha='none')

VPreg_yAB.col.sub = VPAB.alB_reg.mean[VPAB.alB_reg.mean$Cond %in% c('BTW', 'WITHIN'),]

VPreg_yAB.col.hab = VPreg_yAB.col.sub[VPreg_yAB.col.sub$Var %in% c('Habitat'),]
VPreg_yAB.col.clim = VPreg_yAB.col.sub[VPreg_yAB.col.sub$Var %in% c('Climate'),]

xmarg2 <- ggplot(VPreg_yAB.col.hab, aes(x = VP, y=..density.., fill = Cond)) + 
  geom_density(alpha = 0.4)  + xlim(0.05, 0.78) +
  scale_fill_manual(values = c('#BBBBBB','black')) + 
  theme_void() + theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) + 
  theme(legend.position = "none")+ ggtitle('AB model') +
  theme(plot.title = element_text(hjust = 1.2, vjust =  -15, size=14)) 


ymarg2 <- ggplot(VPreg_yAB.col.clim) + 
  geom_density(alpha = 0.4, aes(x = VP, fill = Cond)) + xlim(0.05, 0.78) +
  scale_fill_manual(values = c('#BBBBBB','black')) + 
  theme_void()  + theme(plot.margin = unit(c(0, 0, 0, 0), "cm"))+ 
  theme(legend.position = "none") + 
  coord_flip()

s <- plot_spacer()+theme_void()

# Arranging the plot using cowplot
plot_grid(xmarg, s, xmarg2, s, 
          gPA, ymarg, gAB, ymarg2, ncol = 4, align = "hv", labels = c( ' ', 'a)',' ', 'b)'),
          label_size = 20, #vjust = -0.5,
          rel_widths = c(2, 1), rel_heights = c(1, 2))

library(colorblindr)
cvd_grid(plot_grid(xmarg, s, xmarg2, s, 
                   gPA, ymarg, gAB, ymarg2, ncol = 4, align = "hv", labels = c( ' ', 'a)',' ', 'b)'), 
                   label_size = 20, #vjust = -0.5,
                   rel_widths = c(2, 1), rel_heights = c(1, 2)))


plot_grid(s, s, s, s, 
          gPA, s, gAB, s, ncol = 4, align = "hv", 
          rel_widths = c(2, 1), rel_heights = c(1, 2))

ggarrange(xmarg, s, xmarg2, s, 
         gPA, ymarg, gAB, ymarg2, ncol = 4, nrow=2, widths = c(2,1), heights = c(1,2))
