### TRENDS AND DRIVERS OF SEAGRASSES IN LOWER 
### CHESAPEAKE BAY FROM TRANSECT SURVEYS

# Authors: Lefcheck.Richardson.Orth 
# Last updated: 16 March 2018
# MEPS ms - Warming temperatures alter the relative abundanceand distribution of two
#              co-occurring foundational seagrasses in Chesapeake Bay, USA

library(dplyr)
library(ggplot2)
library(gridExtra)
library(mgcv)
library(plotrix)
library(tidyr)

setwd("R:/Current Projects/Transects/Paul's Files/Master Analysis")
###
# Read in grass dataset
grass <- read.csv("CBP.Zm.Rm.Richardson.csv", header = TRUE)

source("EvaluateSmooths.R")

### PLOTTING TRENDS

# Identify predominantly mesohaline sites
table(grass[which(grass$SALINITY < 15), c("year", "Transect")])

# Generate column for number of quadrats per transect
grass <- grass %>%
  
  group_by(Transect, year) %>%
  
  summarize(Transect.N = length(na.omit(Zm))) %>%
  
  left_join(grass, .)
  
# Summarize data for plotting
grass.summary <- grass %>%
  
  # Remove mesohaline sites
  filter(!Transect %in% c("Dameron_E", "Dameron_S")) %>%

  # Remove restoration sites
  filter(!Transect %in% c("Gwynns", "USCG", "VA_James")) %>%
  
  # Get % cover by transect (no error)
  group_by(year, Transect) %>%
  
  summarize(
    Transect.Area = length(na.omit(Zm)) * quad.size[1],
    Zm = sum(Zm * quad.size, na.rm = T) / Transect.Area,
    Rm = sum(Rm * quad.size, na.rm = T) / Transect.Area,
    depth = mean(depth, na.rm = T) # diff(range(depth, na.rm = T))
  ) 
  
# Take mean & std.error of whole transect % cover
grass.summary2 <- grass.summary %>% 
  
  group_by(year) %>%
  
  summarize(
    Zm.mean = mean(Zm, na.rm = T),
    Zm.se = std.error(Zm, na.rm = T),
    Zm.N = length(na.omit(Zm)),
    Rm.mean = mean(Rm, na.rm = T),
    Rm.se = std.error(Rm, na.rm = T),
    Rm.N = length(na.omit(Rm))
  ) %>%
  
  gather(variable, mean.value, Zm.mean, Rm.mean) %>%
  
  gather(variable2, se.value, Zm.se, Rm.se) %>%
  
  gather(variable3, N, Zm.N, Rm.N) %>%
  
  # Add missing years to generate spaces in plot
  rbind(
    expand.grid(
      year = c(1979, 1991, 2001),
      variable = c("Zm.mean", "Rm.mean"),
      mean.value = NA,
      variable2 = NA,
      se.value = NA,
      variable3 = NA,
      N = NA)
  )

# # Write fn for pooled standard errors
# pooled.se <- function(s, n) {
# 
#   sp <- sqrt((sum((n - 1) * s^2)) / (sum(n) - length(n)))
# 
#   sp * sqrt(sum(1/n))
# 
# }
# 
# # Summarize data for plotting
# grass.summary2 <- grass %>%
# 
#   # Remove mesohaline sites
#   filter(!Transect %in% c("Dameron_E", "Dameron_S")) %>%
# 
#   # Remove restoration sites
#   filter(!Transect %in% c("Gwynns", "USCG", "VA_James")) %>%
# 
#   group_by(year, Transect) %>%
# 
#   summarize(
#     Zm.mean = mean(Zm, na.rm = T),
#     Zm.sd = sd(Zm, na.rm = T),
#     Zm.N = length(Zm),
#     Rm.mean = mean(Rm, na.rm = T),
#     Rm.sd = sd(Rm, na.rm = T),
#     Rm.N = length(Rm)) %>%
# 
#   group_by(year) %>%
# 
#   summarize(
#     Zm.mean = mean(Zm.mean),
#     Zm.se.pooled = pooled.se(Zm.sd, Zm.N),
#     Rm.mean = mean(Rm.mean),
#     Rm.se.pooled = pooled.se(Rm.sd, Rm.N)
#   ) %>%
# 
#   gather(variable, mean.value, Zm.mean, Rm.mean) %>%
# 
#   gather(variable2, se.value, Zm.se.pooled, Rm.se.pooled) %>%
# 
#   # Add missing years to generate spaces in plot
#   rbind(
#     expand.grid(
#       year = c(1979, 1991, 2001),
#       variable = c("Zm.mean", "Rm.mean"),
#       mean.value = NA,
#       variable2 = NA,
#       se.value = NA)
#   )

# Plot results
grass.summary2$variable = as.factor(grass.summary2$variable)

levels(grass.summary2$variable) = c("Ruppia", "Zostera")

grass.summary2[which((grass.summary2$mean.value - grass.summary2$se.value) < 0), "se.value"] <- 0

(pTrend <- ggplot(grass.summary2, 
       aes(x = as.factor(year), y = mean.value, group = variable, col = variable)) +
  geom_line(lwd = 1) +
  geom_errorbar(aes(ymax = mean.value + se.value, ymin = mean.value - se.value),
                width = 0, lwd = 0.8) +
  geom_point(size = 4) +
  geom_vline(xintercept = c(2, 4, 13), lty = 2, col = "grey50") + #c("1979", "1991", "2001")) +
    geom_text(data = subset(grass.summary2, variable == "Ruppia"), 
              aes(label = N, y = -Inf), col = "black", size = 3, vjust = -1) +
  scale_x_discrete(
    labels = c("'78", "", "'91", "", paste0("'", 93:99), "'00", "", paste0("'0", 6:9), paste0("'", 10:16))
    ) +
  ylim(values = c(-2.5, 63)) +
  scale_color_manual(values = c("green2", "forestgreen"), name = "") +
  labs(x = "", y = "Mean % Cover") +
  theme_bw(base_size = 18) +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.background = element_blank(),
    legend.position = c(0.8, 0.9)
  )
)

pdf("Figure 2.pdf",
    width = 8, height = 5)
bquiet = print(pTrend)
dev.off()

# Plot individual trajectories
grass.summary3 = grass.summary %>%
  
  filter(year %in% 2006:2016) %>%
  
  gather(species, cover, Zm, Rm)

grass.summary3$Transect = gsub("\\_", " ", grass.summary3$Transect)

(tTrend = ggplot(grass.summary3, aes(x = year, y = cover, col = species, group = species)) +
  geom_line(lwd = 0.5) +
  geom_point(size = 2) +
  scale_color_manual(values = c("green2", "forestgreen"), name = "") +
  facet_grid(Transect ~ .) + #, scales = "free_y") +
  scale_x_continuous(breaks = c(2006, 2011, 2016)) +
  scale_y_continuous(breaks = c(0, 25, 50)) +
  labs(x = "", y = "Cover") +
  theme_bw(base_size = 18) +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.spacing = unit(0, "lines"),
    axis.text.y = element_text(size = 7),
    strip.background = element_blank(),
    strip.text.y = element_text(angle = 0, hjust = 0, size = 10),
    legend.position = "bottom"
  )
)

pdf("Figure 3.pdf",
    width = 5, height = 12)
bquiet = print(tTrend)
dev.off()

# Look at summary statistics for both species
range(subset(grass.summary2, variable == "Zostera")$mean.value, na.rm = T)

mean(subset(grass.summary2, variable == "Zostera" & year %in% c(1991:2000))$mean.value, na.rm = T)

mean(subset(grass.summary2, variable == "Zostera" & year %in% c(2006:2016))$mean.value, na.rm = T) 

range(subset(grass.summary2, variable == "Ruppia")$mean.value, na.rm = T)

mean(subset(grass.summary2, variable == "Ruppia")$mean.value, na.rm = T)

mean(subset(grass.summary2, variable == "Ruppia" & year %in% c(1991:2000))$mean.value, na.rm = T)

mean(subset(grass.summary2, variable == "Ruppia" & year %in% c(2006:2016))$mean.value, na.rm = T)       

# Look at average cover per transect
mean(subset(grass.summary3, species == "Zm")$cover, na.rm = T)

range(subset(grass.summary3, species == "Zm")$cover, na.rm = T)

mean(subset(grass.summary3, species == "Rm")$cover, na.rm = T)

range(subset(grass.summary3, species == "Rm")$cover, na.rm = T)

# Look at proportion of transects that are dominated by Zostera, Ruppia, or both
# over the course of the survey

grass.summary4 = grass.summary %>% 
  
  filter(year %in% 2006:2016) %>%
  
  group_by(Transect) %>%
  
  summarize(Zm = mean(Zm), Rm = mean(Rm), ratio = Zm/Rm) 

table(cut(grass.summary4$ratio, c(0, 0.5, 1.5, max(grass.summary4$ratio))))

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

### GAM MODELING

# Collapse environmental data
env <- grass %>%
  
  filter(year %in% 2006:2016) %>%
  
  select(year, Transect, CHLA, SALINITY, SECCHI, TN, TP, WTEMP, WTEMP.MAX, X.cord, Y.cord) %>%
  
  group_by(year, Transect) %>%
  
  summarize_each(funs(mean))

# Merge grass summary and environmental data
grass.summary.env <- merge(grass.summary, env)

# Remove years without data
removerows <- apply(grass.summary.env, 1, function(x) all(!is.na(x[4:length(x)])))

grass.summary.env <- grass.summary.env[removerows, ]

# Look at pairwise correlations
cor(grass.summary.env[, (3:13)], use = "complete.obs")

# Run GAM predicting Zm cover
Zm.gam <- gam(Zm ~ 
      s(X.cord, Y.cord, k = 28) +
      s(Rm) + 
      s(CHLA) + 
      s(SALINITY) + 
      s(SECCHI) + 
      s(TN) + 
      s(TP) + 
      s(WTEMP) +
      # s(WTEMP.MAX) +
      s(depth) +
      Transect.Area,
    correlation = corAR1(year ~ 1),
    method = "REML",
    data = grass.summary.env
)
       
summary(Zm.gam)

# plot(Zm.gam)

# Check assumptions
gam.check(Zm.gam)

# Plot gam output
Zm.fitted <- EvaluateSmooths(Zm.gam, select = c(2, 3, 6, 8, 9))

# Plot predictions
(pZm.a = ggplot(subset(Zm.fitted, x.var == "Rm"), aes(x = x.val, y = value)) +
    # geom_hline(yintercept = 3.4e5, alpha = 0) +
    geom_hline(yintercept = 0) +
    geom_ribbon(aes(ymax = value + 2*se, ymin = value - 2*se), fill = "green4", alpha = 0.3) +
    geom_line(col = "green2", lwd = 1.5) +
    # scale_y_continuous(breaks = c(0, 0.2)) +
    annotate("text", x = -Inf, y = Inf, label = "A", hjust = -2, vjust = 1.5, size = 7) +
    labs(x = expression(paste(italic("R. maritima"), "% cover")), 
         y = expression(paste(italic("s"), "(", italic("R. maritima"), ")"))) +
    theme_bw(base_size = 14) +
    theme(
      # plot.margin = unit(c(1, 1, 0, 0), "mm"),
      # panel.background = element_rect(fill = "grey95"),
      panel.grid.major = element_blank(), 
      panel.grid.minor = element_blank(),
      axis.text.y = element_text(angle = 90, hjust = 0.5)
    )
) 

# Plot predictions
(pZm.b = ggplot(subset(Zm.fitted, x.var == "WTEMP"), aes(x = x.val, y = value)) +
    # geom_hline(yintercept = 3.4e5, alpha = 0) +
    geom_hline(yintercept = 0) +
    geom_ribbon(aes(ymax = value + 2*se, ymin = value - 2*se), fill = "firebrick1", alpha = 0.3) +
    geom_line(col = "red", lwd = 1.5) +
    # scale_y_continuous(breaks = c(0, 0.2)) +
    annotate("text", x = -Inf, y = Inf, label = "B", hjust = -2, vjust = 1.5, size = 7) +
    labs(x = expression("Temperature " ( degree*C)), y = expression(paste(italic("s"), "(Temperature)"))) +
    theme_bw(base_size = 14) +
    theme(
      # plot.margin = unit(c(1, 1, 0, 0), "mm"),
      # panel.background = element_rect(fill = "grey95"),
      panel.grid.major = element_blank(), 
      panel.grid.minor = element_blank(),
      axis.text.y = element_text(angle = 90, hjust = 0.5)
    )
) 

# Plot predictions
(pZm.c = ggplot(subset(Zm.fitted, x.var == "TN"), aes(x = x.val, y = value)) +
    # geom_hline(yintercept = 3.4e5, alpha = 0) +
    geom_hline(yintercept = 0) +
    geom_ribbon(aes(ymax = value + 2*se, ymin = value - 2*se), fill = "grey", alpha = 0.3) +
    geom_line(col = "dodgerblue", lwd = 1.5) +
    # scale_y_continuous(breaks = c(0, 0.2)) +
    annotate("text", x = -Inf, y = Inf, label = "C", hjust = -2, vjust = 1.5, size = 7) +
    labs(x = "Total Nitrogen (ug/L)", y = expression(paste(italic("s"), "(Total Nitrogen)"))) +
    theme_bw(base_size = 14) +
    theme(
      # plot.margin = unit(c(1, 1, 0, 0), "mm"),
      # panel.background = element_rect(fill = "grey95"),
      panel.grid.major = element_blank(), 
      panel.grid.minor = element_blank(),
      axis.text.y = element_text(angle = 90, hjust = 0.5)
    )
) 

pdf("Figure 4.pdf",
    width = 8, height = 8)
bquiet = print(
  grid.arrange(pZm.a, pZm.b, pZm.c, nrow = 2)
)
dev.off()

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

# Run GAM predicting Zm cover
Rm.gam <- gam(Rm ~ 
                s(X.cord, Y.cord, k = 28) +
                s(Zm) + 
                s(CHLA) + 
                s(SALINITY) + 
                s(SECCHI) + 
                s(TN) + 
                s(TP) + 
                s(WTEMP) + 
                s(depth) +
                Transect.Area,
              correlation = corAR1(year ~ 1),
              method = "REML",
              data = grass.summary.env
)

summary(Rm.gam)

# plot(Rm.gam)

# Check assumptions
gam.check(Rm.gam)

# Plot gam output
Rm.fitted <- EvaluateSmooths(Rm.gam, select = c(2:4, 8, 9))

# Plot predictions
(pRm.a = ggplot(subset(Rm.fitted, x.var == "Zm"), aes(x = x.val, y = value)) +
    # geom_hline(yintercept = 3.4e5, alpha = 0) +
    geom_hline(yintercept = 0) +
    geom_ribbon(aes(ymax = value + 2*se, ymin = value - 2*se), fill = "forestgreen", alpha = 0.3) +
    geom_line(col = "darkgreen", lwd = 1.5) +
    # scale_y_continuous(breaks = c(0, 0.2)) +
    annotate("text", x = -Inf, y = Inf, label = "A", hjust = -2, vjust = 1.5, size = 7) +
    labs(x = expression(paste(italic("Z. marina"), "% cover")), 
         y = expression(paste(italic("s"), "(", italic("Z. marina"), ")"))) +
    theme_bw(base_size = 14) +
    theme(
      # plot.margin = unit(c(1, 1, 0, 0), "mm"),
      # panel.background = element_rect(fill = "grey95"),
      panel.grid.major = element_blank(), 
      panel.grid.minor = element_blank(),
      axis.text.y = element_text(angle = 90, hjust = 0.5)
    )
) 

# Plot predictions
(pRm.b = ggplot(subset(Rm.fitted, x.var == "CHLA"), aes(x = x.val, y = value)) +
    # geom_hline(yintercept = 3.4e5, alpha = 0) +
    geom_hline(yintercept = 0) +
    geom_ribbon(aes(ymax = value + 2*se, ymin = value - 2*se), fill = "darkorange", alpha = 0.3) +
    geom_line(col = "darkorange3", lwd = 1.5) +
    # scale_y_continuous(breaks = c(0, 0.2)) +
    annotate("text", x = -Inf, y = Inf, label = "B", hjust = -2, vjust = 1.5, size = 7) +
    labs(x = paste("Chl-a"), 
         y = expression(paste(italic("s"), "(Chl-a)"))) +
    theme_bw(base_size = 14) +
    theme(
      # plot.margin = unit(c(1, 1, 0, 0), "mm"),
      # panel.background = element_rect(fill = "grey95"),
      panel.grid.major = element_blank(), 
      panel.grid.minor = element_blank(),
      axis.text.y = element_text(angle = 90, hjust = 0.5)
    )
) 

# Plot predictions
(pRm.c = ggplot(subset(Rm.fitted, x.var == "SALINITY"), aes(x = x.val, y = value)) +
    # geom_hline(yintercept = 3.4e5, alpha = 0) +
    geom_hline(yintercept = 0) +
    geom_ribbon(aes(ymax = value + 2*se, ymin = value - 2*se), fill = "grey", alpha = 0.3) +
    geom_line(col = "black", lwd = 1.5) +
    # scale_y_continuous(breaks = c(0, 0.2)) +
    annotate("text", x = -Inf, y = Inf, label = "C", hjust = -2, vjust = 1.5, size = 7) +
    labs(x = "Salinity", y = expression(paste(italic("s"), "(Salinity)"))) +
    theme_bw(base_size = 14) +
    theme(
      # plot.margin = unit(c(1, 1, 0, 0), "mm"),
      # panel.background = element_rect(fill = "grey95"),
      panel.grid.major = element_blank(), 
      panel.grid.minor = element_blank(),
      axis.text.y = element_text(angle = 90, hjust = 0.5)
    )
) 

pdf("Figure 5.pdf",
    width = 8, height = 8)
bquiet = print(
  grid.arrange(pRm.a, pRm.b, pRm.c, nrow = 2)
)
dev.off()

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

# Plot difference between Zostere and Ruppia against temperature
ZmRmDiff <- grass.summary.env %>%
  
  filter(year > 2007) %>%
  
  group_by(year, Transect, WTEMP) %>%
  
  summarize(diff. = Zm - Rm)

# Plot results
(pDiff <- ggplot(ZmRmDiff, aes(x = WTEMP, y = diff.)) +
  geom_hline(yintercept = 0, lwd = 1, col = "grey50") +
  geom_point(size = 2, shape = 1, col = "grey30") +
  stat_smooth(method = "lm", lwd = 1.4, fill = "forestgreen", col = "darkgreen") +
  labs(x = expression("Temperature " ( degree*C)), 
       y = bquote(italic(Zostera)~~-~italic(Ruppia))) +
  theme_bw(base_size = 18) +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    legend.background = element_blank(),
    legend.position = c(0.8, 0.9)
  )
)

pdf("Figure 6.pdf",
    width = 5, height = 4.5)
bquiet = print(pDiff)
dev.off()

# Significant of LM
summary(lm(diff. ~ WTEMP, ZmRmDiff))

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

# Re-run GAMs ommitting sites with fewer than 3 time points
grass.summary.env2 <- grass.summary.env %>% filter(!Transect %in% c("4_Points_West", "Poquoson_1", "Poquoson_2")) 

Zm.gam2 <- update(Zm.gam, data = grass.summary.env2)

summary(Zm.gam2) # does not make a difference

Rm.gam2 <- update(Rm.gam, data = grass.summary.env2)

summary(Rm.gam2) # does not make a difference