
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)


#......................................................New plotting

setwd("D:/Helsinki/Moth_project/ResultsAug23/RDATA/3.ModelRun")
load('VPycdd-PAclust_1108.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 <- valcondM1H
  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 <- valcondM1H.up
  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 <- valcondM1H.low
  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,M1Ha.vp, along = 3)
  
  M1H.valbw <- aperm(M1Hbw.vp2, c(3,1,2))
  rownames(M1H.valbw) = c('BTW', 'WITHIN', 'ALL')
  
  valbwM1H = as.data.frame(apply(M1H.valbw, 2, c))
  
  valbwM1H$Cond = as.vector(rep(rownames(M1H.valbw), dim(M1Hbw.vp)[3]))
  
  valbw.M1H_long <- valbwM1H
  valbw.M1H_long <- reshape2::melt(valbw.M1H_long)
  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, M1Ha.up, along = 3)
  
  M1H.valbwup <- aperm(M1Hbw.up2, c(3,1,2))
  rownames(M1H.valbwup) = c('BTW', 'WITHIN', 'ALL')
  
  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 <- valbwM1H.up
  
  #### low
  dim(M1Ha.low)
  
  M1H.valbwlow.tmp <- aperm(M1Hbw.low, c(1,3,2))
  dim(M1H.valbwlow.tmp)
  M1Hbw.low2 = abind(M1H.valbwlow.tmp, M1Ha.low, along = 3)
  dim(M1Hbw.low2)
  
  M1H.valbwlow <- aperm(M1Hbw.low2, c(3,1,2))
  rownames(M1H.valbwlow) = c('BTW', 'WITHIN', 'ALL')
  
  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 <- valbwM1H.low
  
  
  # combine
  VP.all1 = rbind(valcond.M1H_long, valbw.M1H_long)
  head(VP.all1)
  
  VP.all1$Cond = factor(VP.all1$Cond, levels=c('Water', 'Forest', 'Wetland - vegetation', 'Heterogeneous', 'BTW', 'WITHIN', 'ALL'))
  

}

VPPA.alB_reg = VP.all1





load('VPycdd-ABclust_1108.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 <- valcondM1H
  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 <- valcondM1H.up
  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 <- valcondM1H.low
  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,M1Ha.vp, along = 3)
  
  M1H.valbw <- aperm(M1Hbw.vp2, c(3,1,2))
  rownames(M1H.valbw) = c('BTW', 'WITHIN', 'ALL')
  
  valbwM1H = as.data.frame(apply(M1H.valbw, 2, c))
  
  valbwM1H$Cond = as.vector(rep(rownames(M1H.valbw), dim(M1Hbw.vp)[3]))
  
  valbw.M1H_long <- valbwM1H
  valbw.M1H_long <- reshape2::melt(valbw.M1H_long)
  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, M1Ha.up, along = 3)
  
  M1H.valbwup <- aperm(M1Hbw.up2, c(3,1,2))
  rownames(M1H.valbwup) = c('BTW', 'WITHIN', 'ALL')
  
  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 <- valbwM1H.up
  
  #### low
  dim(M1Ha.low)
  
  M1H.valbwlow.tmp <- aperm(M1Hbw.low, c(1,3,2))
  dim(M1H.valbwlow.tmp)
  M1Hbw.low2 = abind(M1H.valbwlow.tmp, M1Ha.low, along = 3)
  dim(M1Hbw.low2)
  
  M1H.valbwlow <- aperm(M1Hbw.low2, c(3,1,2))
  rownames(M1H.valbwlow) = c('BTW', 'WITHIN', 'ALL')
  
  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 <- valbwM1H.low
  
  
  # combine
  VP.all1 = rbind(valcond.M1H_long, valbw.M1H_long)
  head(VP.all1)
  
  VP.all1$Cond = factor(VP.all1$Cond, levels=c('Water', 'Forest', 'Wetland - vegetation', 'Heterogeneous', 'BTW', 'WITHIN', 'ALL'))
  
}

VPAB.alB_reg = VP.all1





############################################### Species info

#setwd("~/FilesforR")
setwd("D:/Helsinki/Moth_project/ResultsAug23/RDATA/2.ObsInfo")
Moth <- read.csv("Copy_mothMonitoring.csv", header=T)
head(Moth)
summary(Moth$YearStart)  # 1993-2018 so far

library(maps)
library(mapdata)

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


library(sp)
#library(rgdal)
library(raster)
library(viridis)
# WGS84
CRS("+init=epsg:4326")
CRS("+init=epsg:3067")

MothSites = Moth[, c(5:7)]
Msites= MothSites[!duplicated(MothSites),]

Moth.spdf <- SpatialPointsDataFrame(Msites[,c(3,2)], Msites)
proj4string(Moth.spdf) <-CRS("+proj=utm +zone=35 +ellps=GRS80 +units=m +no_defs")

#plot(Moth.spdf, main="Map of Moths Locations")

newMoth<-spTransform(Moth.spdf, CRS("+init=epsg:4326"))
plot(newMoth, main="Map of Moths Locations")


############################################## VP from model



load('MothclusterFIN.RDATA')
head(Moth_knn)
colnames(Moth_knn)[8] = 'Site'
colnames(Moth_knn)[7] = 'cluster'

Site = as.matrix(newMoth$SiteID)
colnames(Site) = 'Site'
sitebor = merge(Site, Moth_knn, by='Site')

colors = c("cyan3", 'forestgreen',"green3", "brown2", 'yellow', 'darkorange')
colBR = data.frame(Cond=c('Water', 'Forest', 'Wetland - vegetation', 'Heterogeneous', 'BTW', 'WITHIN'), col=colors)

VPPA.alB_reg.mean = VPPA.alB_reg[VPPA.alB_reg$Cond != 'ALL',]
VPAB.alB_reg.mean = VPAB.alB_reg[VPAB.alB_reg$Cond != 'ALL',]

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

# year

####
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)) +
#   geom_point(aes(col=Cond, size=2)) +
#   scale_color_manual(values = c('Heterogeneous'="brown2", 'Forest'="forestgreen", 
#                                 'Water'="cyan3", 'Wetland - vegetation'="green3", 
#                                 'BTW'='yellow', 'WITHIN'='darkorange'))+
#   scale_fill_manual(values = c("brown2", "forestgreen", "darkcyan", "azure3", 'yellow','darkorange'))+
#   theme_pubr(legend = c(0.055,0.78))+
#   ylim(0.07, 0.5) + xlim(0.2, 0.5) +
#   labs(
#     x = "VP habitat",
#     y = "VP climate",
#     fill = "Cluster"
#   ) + geom_abline(intercept = 0, slope = 1) 


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



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