library(bayesplot)  # diagnostic plots of Bayesian models
library(brms)       # Bayesian regression models
library(digest)     # Hashing the data
library(ggbeeswarm) # For bee swarm plot
library(knitr)      # nice tables
library(readODS)    # Reading in open data sheet files
library(tidyverse)  # all things tidy

# Set mcmc options for number of cores to use
options(mc.cores = 15)

# Read in the data
dat <- read_ods(here::here("Data/NSFCSEESEcoData.ods"), "MasterList")

# Authenticate data integrity
dat.hash <- "c371cabbc6c0a8f164e59cc042270438"
if(digest::digest(dat, algo = "md5") != dat.hash) {
  stop("Hashes do not match. Input data may have been changed.")
}

# Make sure the data is in the correct order.
dat <- dat %>% arrange(code)

# Create the local z-score function
local.z <- function(val, valsd) {
  ls.val  <- val[1:13]
  ls.sd   <- valsd[1:13]
  nm.val  <- val[14:26]
  nm.sd   <- valsd[14:26]
  output  <- list()
  pair    <- 1:13
  dif     <- ls.val - nm.val
  z.val   <- dif/sqrt((ls.sd^2 + nm.sd^2)/2)
  
  output$`Pair-level z-scores` <- data.frame(
    Pair = pair,
    zScore = z.val
  )
  
  output$`Overall z-score` <- mean(z.val)
  
  output$`Explanation` <- data.frame(
    zScore = c("Positive: LS > NM",
               "Negative: LS < NM")
  )
  
  return(output)
}

# Create the regional z-score function
regional.z <- function(val) {
  ls.val  <- val[1:13]
  ls.var  <- var(ls.val)
  nm.val  <- val[14:26]
  nm.var  <- var(nm.val)
  output  <- list()
  pair    <- 1:13
  dif     <- ls.val - nm.val
  z.val   <- dif/sqrt((ls.var + nm.var)/2)
  
  output$`Pair-level z-scores` <- data.frame(
    Pair = pair,
    zScore = z.val
  )
  
  output$`Overall z-score` <- mean(z.val)
  
  output$`Explanation` <- data.frame(
    zScore = c("Positive: LS > NM",
               "Negative: LS < NM")
  )
  
  return(output)
}

# Calculate the local z-scores
mus.z <- local.z(dat$Mussels, dat$MusselsSD)
oys.z <- local.z(dat$Oysters, dat$OystersSD)
lit.z <- local.z(dat$Littoraria, dat$LittorariaSD)
bur.z <- local.z(dat$Burrows, dat$BurrowsSD)
alt.z <- local.z(dat$Spartina, dat$SpartinaSD)
org.z <- local.z(dat$OM, dat$OMsd)
car.z <- local.z(dat$C, dat$Csd)
nit.z <- local.z(dat$N, dat$Nsd)
pho.z <- local.z(dat$P, dat$Psd)
lmm.z <- local.z(dat$LMOnly, dat$LMSD)

# Calculate the regional z-scores
fib.z <- regional.z(dat$FishBiomass)
crb.z <- regional.z(dat$Crabs)
shr.z <- regional.z(dat$Shrimp)
fia.z <- regional.z(dat$FishAbu)
jfa.z <- regional.z(dat$JuvenileFish_Abu)
ffa.z <- regional.z(dat$ForageFish_Abu)
fid.z <- regional.z(dat$Fish_diversity)
her.z <- regional.z(dat$HeronUse)
ter.z <- regional.z(dat$Terrapin)

# Create a data frame of the pair-level z-scores for each
# metric
zdat.pair <- data.frame(
  Pair = 1:13,
  Mussels = mus.z$`Pair-level z-scores`$zScore,
  Oysters = oys.z$`Pair-level z-scores`$zScore,
  Littoraria = lit.z$`Pair-level z-scores`$zScore,
  Burrows = bur.z$`Pair-level z-scores`$zScore,
  Spartina = alt.z$`Pair-level z-scores`$zScore,
  OrganicMatter = org.z$`Pair-level z-scores`$zScore,
  Carbon = car.z$`Pair-level z-scores`$zScore,
  Nitrogen = nit.z$`Pair-level z-scores`$zScore,
  Phosphorus = pho.z$`Pair-level z-scores`$zScore,
  FishBio = fib.z$`Pair-level z-scores`$zScore,
  Crabs = crb.z$`Pair-level z-scores`$zScore,
  Shrimp = shr.z$`Pair-level z-scores`$zScore,
  FishAbundance = fia.z$`Pair-level z-scores`$zScore,
  JuvFishAbund = jfa.z$`Pair-level z-scores`$zScore,
  ForageFishAbund = ffa.z$`Pair-level z-scores`$zScore,
  FishDiversity = fid.z$`Pair-level z-scores`$zScore,
  Herons = her.z$`Pair-level z-scores`$zScore,
  Terrapin = ter.z$`Pair-level z-scores`$zScore
)

# Calculate the overall Z score for each pair
zdat.pair$Overall <- rowMeans(zdat.pair[,-1])

# Data frame of the overall scores
zdat.mean <- data.frame(
  Metric = names(zdat.pair)[2:(ncol(zdat.pair)-1)],
  Score = c(
    mus.z$`Overall z-score`,
    oys.z$`Overall z-score`,
    lit.z$`Overall z-score`,
    bur.z$`Overall z-score`,
    alt.z$`Overall z-score`,
    org.z$`Overall z-score`,
    car.z$`Overall z-score`,
    nit.z$`Overall z-score`,
    pho.z$`Overall z-score`,
    fib.z$`Overall z-score`,
    crb.z$`Overall z-score`,
    shr.z$`Overall z-score`,
    fia.z$`Overall z-score`,
    jfa.z$`Overall z-score`,
    ffa.z$`Overall z-score`,
    fid.z$`Overall z-score`,
    her.z$`Overall z-score`,
    ter.z$`Overall z-score`
  )
)

# Give the metrics meaningful names
zdat.mean$Metric <- c(
  "Mussels",
  "Oysters",
  "Periwinkles",
  "Burrows",
  "Cordgrass",
  "Organic Matter",
  "Carbon",
  "Nitrogen",
  "Phosphorus",
  "Fish Biomass",
  "Crab Biomass",
  "Shrimp Biomass",
  "Fish Abundance",
  "Juvenile Fish Abundance",
  "Forage Fish Abundance",
  "Fish Diversity",
  "Heron Use",
  "Terrapin"
)

# What is the mean z-score across all z-scores?
overall.zscore <- mean(zdat.mean$Score)
overall.zscore.sd <- sd(zdat.mean$Score)


# Z-scores vs. age --------------------------------------------------------

# Data frame of overall z-score for each pair with the 
# age of the living shoreline in 2018
sva <- data.frame(
  Overall = zdat.pair$Overall,
  Age = 2018 - dat$Year_Bulit[1:13]
)

# The section below is commented out to speed up manuscript 
# knitting time. Select all commented lines below and press
# CTRL + SHIFT + C to uncomment them for independent
# assessment.

# # Set the priors for the null model
# prior.null <- prior(normal(0,1), class=Intercept) +
#   prior(gamma(1,1), class = sigma)
# # Specify the null model
# m0 <- brm(
#   Overall ~ 1,
#   data = sva,
#   family = gaussian(),
#   prior = prior.null,
#   seed = 2276,
#   warmup = 5e3,
#   iter = 5e4,
#   control = list(adapt_delta = 0.95),
#   save_pars = save_pars(all = TRUE)
# )
# # Save null model
# saveRDS(
#   object = m0,
#   file = here::here("Output/Models/AgeNull.rds")
# )

# Read in the null model
m0 <- readRDS(here::here("Output/Models/AgeNull.rds"))

# #Set the priors for the age model
# prior1 <- prior(normal(0,1), nlpar = "b0") +
#   prior(normal(0.25, 0.75), nlpar = "b1") +
#   prior(gamma(1,1), class = sigma)
# m1 <- brm(
#   bf(Overall ~ b0 + b1*log(Age),
#      b0 + b1 ~ 1,
#      nl = TRUE),
#   data = sva,
#   family = gaussian(),
#   prior = prior1,
#   seed = 2276,
#   warmup = 5e3,
#   iter = 5e4,
#   control = list(adapt_delta = 0.95),
#   save_pars = save_pars(all = TRUE)
# )
# # Save age model
# saveRDS(
#   object = m1,
#   file = here::here("Output/Models/AgeModel.rds")
# )

# Read in the age model
m1 <- readRDS(here::here("Output/Models/AgeModel.rds"))

# Summary of m1
m1.sum <- round(as.data.frame(summary(m1)$fixed),2)

# Check the posteriors
# mvar.m1 <- pp_check(m1, type = "stat_2d")

# Compare m1 to the null
# loo.m1 <- loo(m1, m0, moment_match = TRUE)
# saveRDS(
#   object = loo.m1,
#   file = here::here("Output/Models/LOO.rds")
# )

# Read in the LOO comparison
loo.m1 <- readRDS(here::here("Output/Models/LOO.rds"))

# Get the posterior plot for the effect of age
age.post <- bayesplot::mcmc_areas(
  x = m1,
  pars = "b_b1_Intercept",
  prob = 0.95
) +
  scale_y_discrete(labels = "Age")

# Calculate the proportion of the distribution for the age
# coefficient that overlaps zero
post <- posterior_samples(
  x = m1,
  pars = "b_b1_Intercept"
)
post.zero <- length(post[post <= 0])/nrow(post)

# Plot the mean Z-score against living shoreline age.
SVAplot <- sva %>% ggplot(aes(x = Age, y = Overall)) +
  geom_point(size = 3) +
  theme_classic() + 
  ylab(expression(paste("Mean ",italic("Z"),"-score"))) +
  xlab("Living Shoreline Age (c. 2018)") +
  theme(
    axis.text = element_text(size = 12, color = "black"),
    axis.title = element_text(size = 12)
  )

# SVAplot
# ggsave("NSFCSEESEcoSynthesis/Figure4.pdf",
#        SVAplot,
#        device = "pdf",
#        width = 4,
#        height = 3,
#        units = "in",
#        dpi = 600)



# Without the Sill --------------------------------------------------------

zdat.pair2 <- data.frame(
  Pair = 1:13,
  Mussels = lmm.z$`Pair-level z-scores`$zScore,
  Oysters = oys.z$`Pair-level z-scores`$zScore,
  Littoraria = lit.z$`Pair-level z-scores`$zScore,
  Burrows = bur.z$`Pair-level z-scores`$zScore,
  Spartina = alt.z$`Pair-level z-scores`$zScore,
  OrganicMatter = org.z$`Pair-level z-scores`$zScore,
  Carbon = car.z$`Pair-level z-scores`$zScore,
  Nitrogen = nit.z$`Pair-level z-scores`$zScore,
  Phosphorus = pho.z$`Pair-level z-scores`$zScore,
  FishBio = fib.z$`Pair-level z-scores`$zScore,
  Crabs = crb.z$`Pair-level z-scores`$zScore,
  Shrimp = shr.z$`Pair-level z-scores`$zScore,
  FishAbundance = fia.z$`Pair-level z-scores`$zScore,
  JuvFishAbund = jfa.z$`Pair-level z-scores`$zScore,
  ForageFishAbund = ffa.z$`Pair-level z-scores`$zScore,
  FishDiversity = fid.z$`Pair-level z-scores`$zScore,
  Herons = her.z$`Pair-level z-scores`$zScore,
  Terrapin = ter.z$`Pair-level z-scores`$zScore
)

# Calculate the overall Z score for each pair
zdat.pair2$Overall <- rowMeans(zdat.pair2[,-1])

# Data frame of the overall scores
zdat.mean2 <- data.frame(
  Metric = names(zdat.pair2)[2:(ncol(zdat.pair2)-1)],
  Score = c(
    lmm.z$`Overall z-score`,
    oys.z$`Overall z-score`,
    lit.z$`Overall z-score`,
    bur.z$`Overall z-score`,
    alt.z$`Overall z-score`,
    org.z$`Overall z-score`,
    car.z$`Overall z-score`,
    nit.z$`Overall z-score`,
    pho.z$`Overall z-score`,
    fib.z$`Overall z-score`,
    crb.z$`Overall z-score`,
    shr.z$`Overall z-score`,
    fia.z$`Overall z-score`,
    jfa.z$`Overall z-score`,
    ffa.z$`Overall z-score`,
    fid.z$`Overall z-score`,
    her.z$`Overall z-score`,
    ter.z$`Overall z-score`
  )
)

# What is the mean z-score across all z-scores?
overall.zscore2 <- mean(zdat.mean2$Score)
overall.zscore.sd2 <- sd(zdat.mean2$Score)

# Age vs. Score
sva2 <- data.frame(
  Overall = zdat.pair2$Overall,
  Age = 2018 - dat$Year_Bulit[1:13]
)

# # Set the priors for the null model
# nosill.prior.null <- prior(normal(0,1), class=Intercept) +
#   prior(gamma(1,1), class = sigma)
# # Specify the null model
# no.sill.m0 <- brm(
#   Overall ~ 1,
#   data = sva2,
#   family = gaussian(),
#   prior = nosill.prior.null,
#   seed = 2276,
#   warmup = 5e3,
#   iter = 5e4,
#   control = list(adapt_delta = 0.95),
#   save_pars = save_pars(all = TRUE)
# )
# # Save null model
# saveRDS(
#   object = no.sill.m0,
#   file = here::here("Output/Models/NoSillAgeNull.rds")
# )

# Read in the null model
no.sill.m0 <- readRDS(here::here("Output/Models/NoSillAgeNull.rds"))

# #Set the priors for the age model
# no.sill.prior1 <- prior(normal(0,1), nlpar = "b0") +
#   prior(normal(0.25, 0.75), nlpar = "b1") +
#   prior(gamma(1,1), class = sigma)
# no.sill.m1 <- brm(
#   bf(Overall ~ b0 + b1*log(Age),
#      b0 + b1 ~ 1,
#      nl = TRUE),
#   data = sva2,
#   family = gaussian(),
#   prior = no.sill.prior1,
#   seed = 2276,
#   warmup = 5e3,
#   iter = 5e4,
#   control = list(adapt_delta = 0.95),
#   save_pars = save_pars(all = TRUE)
# )
# # Save age model
# saveRDS(
#   object = no.sill.m1,
#   file = here::here("Output/Models/NoSillAgeModel.rds")
# )

# Read in the age model
no.sill.m1 <- readRDS(here::here("Output/Models/NoSillAgeModel.rds"))

# Summary of m1
no.sill.m1.sum <- round(as.data.frame(summary(no.sill.m1)$fixed),2)

# Check the posteriors
# no.sill.mvar.m1 <- pp_check(no.sill.m1, type = "stat_2d")

# # Compare m1 to the null
# no.sill.loo.m1 <- loo(no.sill.m1, no.sill.m0, moment_match = TRUE)
# saveRDS(
#   object = no.sill.loo.m1,
#   file = here::here("Output/Models/NoSillLOO.rds")
# )

# Read in the LOO comparison
no.sill.loo.m1 <- readRDS(here::here("Output/Models/NoSillLOO.rds"))

# Get the posterior plot for the effect of age
no.sill.age.post <- bayesplot::mcmc_areas(
  x = no.sill.m1,
  pars = "b_b1_Intercept",
  prob = 0.95
) +
  scale_y_discrete(labels = "Age")

# Calculate the proportion of the distribution for the age
# coefficient that overlaps zero
no.sill.post <- posterior_samples(
  x = no.sill.m1,
  pars = "b_b1_Intercept"
)
no.sill.post.zero <- length(no.sill.post[no.sill.post <= 0])/nrow(no.sill.post)

# Bee swarm plots ---------------------------------------------------------

vars <- c("Type","Mussels","Oysters","Littoraria","Burrows","Spartina",
          "FishBiomass","Crabs","Shrimp","FishAbu","JuvenileFish_Abu",
          "ForageFish_Abu","Fish_diversity","HeronUse","OM","C","N",
          "P","Terrapin", "LMOnly")

plot.dat <- dat %>%
  select(all_of(vars))

# plot.dat.long <- pivot_longer(plot.dat,
#                               cols = 2:19,
#                               names_to = "Metric",
#                               values_to = "Value")

mus.plot <- plot.dat %>% 
  ggplot(aes(x = Type,
             y = Mussels)) +
  geom_quasirandom(width = 0.15,
                   size = 2,
                   shape = 21) +
  theme_classic() +
  xlab("Type") +
  ylab(bquote('Mussel Density '~('# m'^-2))) +
  theme(
    axis.text = element_text(size = 8, color = "black"),
    axis.title = element_text(size = 8, color = "black")
  )
nsm.plot <- plot.dat %>% 
  ggplot(aes(x = Type,
             y = LMOnly)) +
  geom_quasirandom(width = 0.15,
                   size = 2,
                   shape = 21) +
  theme_classic() +
  xlab("Type") +
  ylab(bquote('LM Mussel Density '~('# m'^-2))) +
  theme(
    axis.text = element_text(size = 8, color = "black"),
    axis.title = element_text(size = 8, color = "black")
  )
oys.plot <- plot.dat %>% 
  ggplot(aes(x = Type,
             y = Oysters)) +
  geom_quasirandom(width = 0.15,
                   size = 2,
                   shape = 21) +
  theme_classic() +
  xlab("Type") +
  ylab(bquote('Oyster Density '~('# m'^-2))) +
  theme(
    axis.text = element_text(size = 8, color = "black"),
    axis.title = element_text(size = 8, color = "black")
  )
lit.plot <- plot.dat %>% 
  ggplot(aes(x = Type,
             y = Littoraria)) +
  geom_quasirandom(width = 0.15,
                   size = 2,
                   shape = 21) +
  theme_classic() +
  xlab("Type") +
  ylab(bquote('Periwinkle Density '~('# m'^-2))) +
  theme(
    axis.text = element_text(size = 8, color = "black"),
    axis.title = element_text(size = 8, color = "black")
  )
bur.plot <- plot.dat %>% 
  ggplot(aes(x = Type,
             y = Burrows)) +
  geom_quasirandom(width = 0.15,
                   size = 2,
                   shape = 21) +
  theme_classic() +
  xlab("Type") +
  ylab(bquote('Burrow Density '~('burrows m'^-2))) +
  theme(
    axis.text = element_text(size = 8, color = "black"),
    axis.title = element_text(size = 8, color = "black")
  )
bla.plot <- ggplot() + theme_void()
car.plot <- plot.dat %>% 
  ggplot(aes(x = Type,
             y = C)) +
  geom_quasirandom(width = 0.15,
                   size = 2,
                   shape = 21) +
  theme_classic() +
  xlab("Type") +
  ylab(bquote('Carbon '~('% dry weight'))) +
  theme(
    axis.text = element_text(size = 8, color = "black"),
    axis.title = element_text(size = 8, color = "black")
  )
nit.plot <- plot.dat %>% 
  ggplot(aes(x = Type,
             y = N)) +
  geom_quasirandom(width = 0.15,
                   size = 2,
                   shape = 21) +
  theme_classic() +
  xlab("Type") +
  ylab(bquote('Nitrogen '~('% dry weight'))) +
  theme(
    axis.text = element_text(size = 8, color = "black"),
    axis.title = element_text(size = 8, color = "black")
  )
pho.plot <- plot.dat %>% 
  ggplot(aes(x = Type,
             y = P)) +
  geom_quasirandom(width = 0.15,
                   size = 2,
                   shape = 21) +
  theme_classic() +
  xlab("Type") +
  ylab(bquote('Phosphorus '~('% dry weight'))) +
  theme(
    axis.text = element_text(size = 8, color = "black"),
    axis.title = element_text(size = 8, color = "black")
  )
org.plot <- plot.dat %>% 
  ggplot(aes(x = Type,
             y = OM)) +
  geom_quasirandom(width = 0.15,
                   size = 2,
                   shape = 21) +
  theme_classic() +
  xlab("Type") +
  ylab(bquote('Organic Matter '~('% AFDW'))) +
  theme(
    axis.text = element_text(size = 8, color = "black"),
    axis.title = element_text(size = 8, color = "black")
  )
her.plot <- plot.dat %>% 
  ggplot(aes(x = Type,
             y = HeronUse)) +
  geom_quasirandom(width = 0.15,
                   size = 2,
                   shape = 21) +
  theme_classic() +
  xlab("Type") +
  ylab(bquote('Heron Use '~('minutes survey'^-1))) +
  theme(
    axis.text = element_text(size = 8, color = "black"),
    axis.title = element_text(size = 8, color = "black")
  )
ter.plot <- plot.dat %>% 
  ggplot(aes(x = Type,
             y = Terrapin)) +
  geom_quasirandom(width = 0.15,
                   size = 2,
                   shape = 21) +
  theme_classic() +
  xlab("Type") +
  ylab(bquote('Terrapin Density '~('# site'^-1))) +
  theme(
    axis.text = element_text(size = 8, color = "black"),
    axis.title = element_text(size = 8, color = "black")
  )
alt.plot <- plot.dat %>% 
  ggplot(aes(x = Type,
             y = Spartina)) +
  geom_quasirandom(width = 0.15,
                   size = 2,
                   shape = 21) +
  theme_classic() +
  xlab("Type") +
  ylab(bquote('Cordgrass Density '~('stems m'^-2))) +
  theme(
    axis.text = element_text(size = 8, color = "black"),
    axis.title = element_text(size = 8, color = "black")
  )
crb.plot <- plot.dat %>% 
  ggplot(aes(x = Type,
             y = Crabs)) +
  geom_quasirandom(width = 0.15,
                   size = 2,
                   shape = 21) +
  theme_classic() +
  xlab("Type") +
  ylab(bquote('Blue Crab Biomass '~('g total'))) +
  theme(
    axis.text = element_text(size = 8, color = "black"),
    axis.title = element_text(size = 8, color = "black")
  )
shr.plot <- plot.dat %>% 
  ggplot(aes(x = Type,
             y = Shrimp)) +
  geom_quasirandom(width = 0.15,
                   size = 2,
                   shape = 21) +
  theme_classic() +
  xlab("Type") +
  ylab(bquote('Shrimp Biomass '~('g total'))) +
  theme(
    axis.text = element_text(size = 8, color = "black"),
    axis.title = element_text(size = 8, color = "black")
  )
ffa.plot <- plot.dat %>% 
  ggplot(aes(x = Type,
             y = ForageFish_Abu)) +
  geom_quasirandom(width = 0.15,
                   size = 2,
                   shape = 21) +
  theme_classic() +
  xlab("Type") +
  ylab(bquote('Forage Fish Abundance '~('# y'^-1))) +
  theme(
    axis.text = element_text(size = 8, color = "black"),
    axis.title = element_text(size = 8, color = "black")
  )
fia.plot <- plot.dat %>% 
  ggplot(aes(x = Type,
             y = FishAbu)) +
  geom_quasirandom(width = 0.15,
                   size = 2,
                   shape = 21) +
  theme_classic() +
  xlab("Type") +
  ylab(bquote('Fish Abundance '~('# total'))) +
  theme(
    axis.text = element_text(size = 8, color = "black"),
    axis.title = element_text(size = 8, color = "black")
  )
fib.plot <- plot.dat %>% 
  ggplot(aes(x = Type,
             y = FishBiomass)) +
  geom_quasirandom(width = 0.15,
                   size = 2,
                   shape = 21) +
  theme_classic() +
  xlab("Type") +
  ylab(bquote('Fish Biomass '~('g total'))) +
  theme(
    axis.text = element_text(size = 8, color = "black"),
    axis.title = element_text(size = 8, color = "black")
  )
fid.plot <- plot.dat %>% 
  ggplot(aes(x = Type,
             y = Fish_diversity)) +
  geom_quasirandom(width = 0.15,
                   size = 2,
                   shape = 21) +
  theme_classic() +
  xlab("Type") +
  ylab(bquote('Fish Diversity '~('TD y'^-1))) +
  theme(
    axis.text = element_text(size = 8, color = "black"),
    axis.title = element_text(size = 8, color = "black")
  )
jfa.plot <- plot.dat %>% 
  ggplot(aes(x = Type,
             y = JuvenileFish_Abu)) +
  geom_quasirandom(width = 0.15,
                   size = 2,
                   shape = 21) +
  theme_classic() +
  xlab("Type") +
  ylab(bquote('Juvenile Fish Abundance '~('# y'^-1))) +
  theme(
    axis.text = element_text(size = 8, color = "black"),
    axis.title = element_text(size = 8, color = "black")
  )

allPlots <- ggpubr::ggarrange(car.plot, nit.plot, pho.plot, org.plot,
                              crb.plot, shr.plot, fib.plot, fia.plot,
                              ffa.plot, jfa.plot, fid.plot, alt.plot,
                              mus.plot, oys.plot, lit.plot, bur.plot,
                              nsm.plot, her.plot, ter.plot,
                              ncol = 4, nrow = 5)
# ggsave(
#   plot = allPlots,
#   filename = here::here("Figures/Figure2.pdf"),
#   width = 6.5,
#   height = 8.5,
#   units = "in",
#   dpi = 600,
#   device = "pdf"
# )


# Z-score plot ------------------------------------------------------------

# Create a data frame of the pair-level z-scores for each
# metric including both mussel metrics
zdat.pair3 <- data.frame(
  Pair = 1:13,
  Mussels = mus.z$`Pair-level z-scores`$zScore,
  Mussels2 = lmm.z$`Pair-level z-scores`$zScore,
  Oysters = oys.z$`Pair-level z-scores`$zScore,
  Littoraria = lit.z$`Pair-level z-scores`$zScore,
  Burrows = bur.z$`Pair-level z-scores`$zScore,
  Spartina = alt.z$`Pair-level z-scores`$zScore,
  OrganicMatter = org.z$`Pair-level z-scores`$zScore,
  Carbon = car.z$`Pair-level z-scores`$zScore,
  Nitrogen = nit.z$`Pair-level z-scores`$zScore,
  Phosphorus = pho.z$`Pair-level z-scores`$zScore,
  FishBio = fib.z$`Pair-level z-scores`$zScore,
  Crabs = crb.z$`Pair-level z-scores`$zScore,
  Shrimp = shr.z$`Pair-level z-scores`$zScore,
  FishAbundance = fia.z$`Pair-level z-scores`$zScore,
  JuvFishAbund = jfa.z$`Pair-level z-scores`$zScore,
  ForageFishAbund = ffa.z$`Pair-level z-scores`$zScore,
  FishDiversity = fid.z$`Pair-level z-scores`$zScore,
  Herons = her.z$`Pair-level z-scores`$zScore,
  Terrapin = ter.z$`Pair-level z-scores`$zScore
)

# Calculate the overall Z score for each pair
zdat.pair3$Overall <- rowMeans(zdat.pair3[,-c(1,3)]) # Exclude Mussels2

# # Save the output
# data.table::fwrite(
#   x = zdat.pair3,
#   file = here::here("Output/Data/ZScoresByPair.csv"))

# Data frame of the overall scores
zdat.mean3 <- data.frame(
  Metric = names(zdat.pair3)[2:(ncol(zdat.pair3)-1)],
  Score = c(
    mus.z$`Overall z-score`,
    lmm.z$`Overall z-score`,
    oys.z$`Overall z-score`,
    lit.z$`Overall z-score`,
    bur.z$`Overall z-score`,
    alt.z$`Overall z-score`,
    org.z$`Overall z-score`,
    car.z$`Overall z-score`,
    nit.z$`Overall z-score`,
    pho.z$`Overall z-score`,
    fib.z$`Overall z-score`,
    crb.z$`Overall z-score`,
    shr.z$`Overall z-score`,
    fia.z$`Overall z-score`,
    jfa.z$`Overall z-score`,
    ffa.z$`Overall z-score`,
    fid.z$`Overall z-score`,
    her.z$`Overall z-score`,
    ter.z$`Overall z-score`
  )
)

# Give the metrics meaningful names
zdat.mean3$Metric <- c(
  "Mussels",
  "Mussels (LM Only)",
  "Oysters",
  "Periwinkles",
  "Burrows",
  "Cordgrass",
  "Organic Matter",
  "Carbon",
  "Nitrogen",
  "Phosphorus",
  "Fish Biomass",
  "Crab Biomass",
  "Shrimp Biomass",
  "Fish Abundance",
  "Juvenile Fish Abundance",
  "Forage Fish Abundance",
  "Fish Diversity",
  "Heron Use",
  "Terrapin"
)

# Pivot longer
zlong <- zdat.pair3 %>% 
  select(-21) %>% # remove the overall score
  pivot_longer(cols = 2:20,
               values_to = "Score",
               names_to = "Metric")
zlong$Metric <- factor(zlong$Metric, levels = unique(zlong$Metric))

# Bee swarm plot
z.bee <- zlong %>% 
  ggplot(aes(x = Metric,
             y = Score)) +
  geom_hline(yintercept = 0, linetype = 2) +
  # geom_vline(xintercept = 1:18, linetype = 3, size = 0.2) +
  geom_quasirandom(width = 0.15,
                   size = 2,
                   shape = 21) +
  theme_classic() +
  xlab("Metric") +
  ylab("Z-score") +
  scale_x_discrete(labels = rev(zdat.mean3$Metric),
                   limits = rev) +
  theme(
    axis.text = element_text(size = 12, color = "black"),
    axis.title = element_text(size = 12, color = "black")
  ) +
  coord_flip()
# ggsave(
#   filename = "Figures/Figure3.pdf",
#   plot = z.bee,
#   device = "pdf",
#   height = 6,
#   width = 5.5,
#   units = "in",
#   dpi = 600
# )


# Save the session info ---------------------------------------------------

sink(here::here("Output/SessionInfo.txt"))
devtools::session_info()
sink()