fishdrought_resubmittal.Rmd

---
title: "Fish multi-yeardrought"
author: "Rosemary Hartman"
date: "`r Sys.Date()`"
output:
  pdf_document: default
  html_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, warning = FALSE)
library(tidyverse)
library(readxl)
library(viridis)
library(lubridate)
library(emmeans)
library(visreg)
library(ggsignif)
library(DroughtData)
library(zoo)
library("ggbeeswarm")
library(broom)
 library(gridExtra)
```

## Fish - multiple drought years

I categorized drought years into year '0' (wet years), '1' (single dry year or first year of drought), '2' (Second year of drought) and '3+' (third or more year of drought. )

```{r echo = F}
load("FishDrough.RData")
#create a "long term drought" index that is the sum of the previous three year's CVI
yrs = read_csv("data/yearassignments.csv")  %>%
  select(Year, DroughtYear2)

FishDrought = left_join(FishDrought, yrs) %>%
  mutate(DroughtYear3 = as.ordered(DroughtYear)) %>%
  group_by(MetricL) %>%
  arrange(Year) %>%
  mutate(LagValue = lag(Value))



ggplot(FishDrought, aes(x = DroughtYear3, y = Value, fill = DroughtYear3)) + 
  geom_boxplot() +
  geom_quasirandom()+
  facet_wrap(MetricL~., scales = "free_y")+ theme_bw()+
  scale_fill_brewer(palette = "Dark2", name = NULL, labels = c("Wet years", "First/only dry year",
                                                                            "Second dry year", "Thrid or more dry year"))+
  theme(legend.position = "bottom")


ggplot(FishDrought, aes(x = Year, y = Value))+ geom_point()+ geom_line()+
  facet_wrap(MetricL~., scales = "free_y")+ theme_bw()

ggplot(FishDrought, aes(x = DroughtYear2, y = Value, color = MetricL)) + 
  geom_point()+
  geom_smooth(method = "lm")

ggplot(FishDrought, aes(x = Index, y = Value, color = MetricL)) + 
  geom_point()+
  geom_smooth(method = "lm")



```

Then I ran statistics. For POD species i included both "year of drought" and "year" to account for changes over time.

### First, striped bass

```{r}
#add a lag term to account for previous year's population
#Or matt suggests N-4 because most striped bass mature at age 4. 
#But N-4 didn't work so I cut it.


SB = filter(FishDrought, Metric == "logSB") 

#Drought year as an ordered factor
FD1 = lm(Value ~ DroughtYear3 + Year + LagValue, data = SB)
acf(residuals(FD1))
summary(FD1)
plot(FD1)
emmeans(FD1, pairwise ~ DroughtYear3)
visreg(FD1)
```


```{r}
#drought year as a non-ordered factor
FD1.1 = lm(Value ~ DroughtYear + Year+ LagValue, data = SB)

summary(FD1.1)
plot(FD1.1)
emmeans(FD1.1, pairwise ~ DroughtYear)
visreg(FD1.1)
```


```{r}
#drought year as a continuous variable
FD1.2 = lm(Value ~ DroughtYear2 + Year+ LagValue, data = SB)
summary(FD1.2)
plot(FD1.2)
visreg(FD1.2)


library(AICcmodavg)
models = list(FD1, FD1.1, FD1.2)
modnames = c("ordered factor", "factor", "continuous")
aictab(models, modnames)
aictab(models)

```
So, droughts are bad, but long droughts aren't much worse than short ones after accounting for long-term population trajectory

OK, visually, what is going on? 

```{r}
ggplot(SB, aes(x = Year, y = Value, fill = DroughtYear))+ geom_col()


```




### Now longfin smelt

```{r}
LFS = filter(FishDrought, Metric == "logLFS")

LFS1 = lm(Value ~ DroughtYear + Year + LagValue, data = LFS)
acf(residuals(LFS1))
summary(LFS1)
plot(LFS1)
emmeans(LFS1, pairwise ~ DroughtYear)
visreg(LFS1)

```

Again, droughts are bad, but long droughts aren't much worse than short ones after accounting for long-term population trajectory

```{r}

LFS1.1 = lm(Value ~ DroughtYear + Year+LagValue, data = LFS)
summary(LFS1.1)
plot(LFS1.1)
emmeans(LFS1.1, pairwise ~ DroughtYear)
visreg(LFS1.1)

LFS1.2 = lm(Value ~ DroughtYear2 + Year+LagValue, data = LFS)
summary(LFS1.2)
plot(LFS1.2)
visreg(LFS1.2)


models = list(LFS1, LFS1.1, LFS1.2)
modnames = c("ordered factor", "factor", "continuous")
aictab(models, modnames)


```

### Delta Smelt

```{r}
DS = filter(FishDrought, Metric == "logDS") 

DS1 = lm(Value ~ DroughtYear + Year + LagValue, data = DS)
acf(residuals(DS1))

summary(DS1)
plot(DS1)
visreg(DS1)
```


```{r}
DS1.1 = lm(Value ~ DroughtYear + Year+LagValue, data = DS)
acf(residuals(DS1.1))
summary(DS1.1)
plot(DS1.1)
visreg(DS1.1)
```


```{r}
DS1.2 = lm(Value ~ DroughtYear2 + Year+ LagValue, data = DS)
summary(DS1.2)
plot(DS1.2)
visreg(DS1.2)

models = list(DS1, DS1.1, DS1.2)

aictab(models, modnames)

```
Delta Smelt just don't want to behave. 

### American Shad

```{r}
shad = filter(FishDrought, Metric == "logShad") 

shad1 = lm(Value ~ DroughtYear + Year+ LagValue, data = shad)
acf(residuals(shad1))
summary(shad1)
emmeans(shad1, pairwise ~ DroughtYear)
plot(shad1)
visreg(shad1)
```


```{r}
shad1.1 = lm(Value ~ DroughtYear + Year+LagValue, data = shad)
summary(shad1.1)
acf(residuals(shad1.1))
emmeans(shad1.1, pairwise ~ DroughtYear)
plot(shad1.1)
visreg(shad1.1)
```


```{r}
shad1.2 = lm(Value ~ DroughtYear2 + Year+LagValue, data = shad)
summary(shad1.2)

plot(shad1.2)
visreg(shad1.2)


```
It's all about dry versus wet, not multi-year droughts

### Now Salmon CRR


I wasn't quite sure if I needed Year as a predictor here, so I did a quick plot

```{r}
salmon = read_csv("data/droughtsalmon.csv")

ggplot(salmon, aes(x = DroughtYearSalmon, y = Value))+ geom_boxplot() +
  facet_wrap(~Metric)
ggplot(salmon, aes(x = DroughtYear, y = Value, fill = DroughtYearSalmon))+ geom_boxplot() +
  facet_wrap(~Metric)

ggplot(salmon, aes(x = Year, y = Value))+ geom_line() +
  facet_wrap(~Metric)

#what about the lag term?


# CRRs = filter(FishDrought, Metric == "Salmon CRR") 
# ggplot(CRRs, aes(x=Year, y = Value))+ geom_point()+ geom_smooth()
# ggplot(CRRs, aes(x=Value, y = LagValue))+ geom_point()+ geom_smooth()

```

No real trend over time, so I'll skip that. 


I need to check with Evan. I'm not sure the lag term makes sense ecologically here. 
``

```{r}
#Spring run
SR = filter(salmon, Metric == "CV SR CRR")
SRcrr = lm(Value ~ DroughtYear, data = SR)
#SRcrr = lm(Value ~ DroughtYearSalmon, data = SR)
summary(SRcrr)
plot(SRcrr)
acf(residuals(SRcrr))
emmeans(SRcrr, pairwise ~ DroughtYear)
visreg(SRcrr)

#WR run
WR = filter(salmon, Metric == "WR CRR")
WRcrr = lm(Value ~ DroughtYear, data = WR)
summary(WRcrr)
plot(WRcrr)
acf(residuals(WRcrr))
emmeans(WRcrr, pairwise ~ DroughtYear)
visreg(WRcrr)
#fall run
FR = filter(salmon, Metric == "CV FR CRR")
FRcrr = lm(Value ~ DroughtYear, data = FR)
summary(FRcrr)
plot(FRcrr)
acf(residuals(FRcrr))
emmeans(FRcrr, pairwise ~ DroughtYear)
visreg(FRcrr)

```


THree or more years of drought is REALLY bad for salmon. This means they are Outmigrating on a droughht year and returning on a drought year. 

Is it more appropriate to give CRR's during outmigraiton year or regurn year? or both?

```{r}



#does it make sense to have the lag of the cohort replacement rate in for salmon? Probably not.

ggplot(CRRs, aes(Value, LagValue))+ geom_point( aes(color = Drought))+ geom_smooth(method = "lm")

```



### Threadfin

```{r}

TFshad = filter(FishDrought, Metric == "logTFS")
TFshad1 = lm(Value ~ DroughtYear + Year + LagValue, data = TFshad)
summary(TFshad1)
emmeans(TFshad1, pairwise ~ DroughtYear)
plot(TFshad1)
visreg(TFshad1)


```



```{r Fig3, warn = FALSE }
#final plot for paper - 

#data frame of pairwise comparison results
# FishDrought = droplevels(FishDrought)
# pairs = data.frame(MetricL = rep(levels(FishDrought$MetricL), each = 4), 
#                    DroughtYear = rep(c("0", "1", "2", "3+"), 5),
#                    Group = c("A", "B", "B", "B", "A", "A", "A", "A", "A", 
#                              "B", "B", "B", "A", "B", "B", "B",
#                              "A", "A", "A", "B"),
#                    Y = c(rep(10, 4), rep(8, 4), rep(12, 4), rep(10, 4), rep(4, 4)))

ggplot(FishDrought, aes(x = as.factor(DroughtYear), y = Value, fill = as.factor(DroughtYear))) + 
  geom_boxplot() +
  geom_quasirandom()+
 # geom_signif(comparisons = list(c("0", "1"), c("0", "2"), c("0", "3+"), c("1", "2"), 
#                                  c("1", "3+"), c("2, 3+")), position = c())+
   theme_bw()+
 # geom_text(data = pairs, aes(y = Y, label = Group))+
  xlab(NULL)+ ylab(NULL)+
  scale_fill_brewer(palette = "Dark2", name = NULL, labels = c("Wet years", "First/only dry year",
                                                                            "Second dry year", "Thrid or more dry year"))+

  theme(legend.position = "bottom")


  
LF = ggplot(filter(FishDrought, Metric == "logLFS"),  aes(x = as.factor(DroughtYear), y = Value, fill = as.factor(DroughtYear))) + 
  geom_boxplot() +
  geom_quasirandom()+
  geom_signif(comparisons = list(c("0", "1"), c("0", "2"), c("0", "3+"), 
                                 c("1", "2"), c("1", "3+"), c("2", "3+")), 
              y_position = c(10,11,12,8, 9, 7), test = t.test, 
              map_signif_level=T)+
  facet_wrap(MetricL~., scales = "free_y")+ theme_bw()+
 # geom_text(data = pairs, aes(y = Y, label = Group))+
  xlab("Drought Length (Years)")+ ylab("log FMWT index")+ ylim(0, 13)+
  scale_fill_brewer(palette = "Dark2", name = NULL, guide = NULL)
LF
```


```{r Fig3a, warn = FALSE }
sh = ggplot(filter(FishDrought, Metric == "logShad"),  aes(x = as.factor(DroughtYear), y = Value, fill = as.factor(DroughtYear))) + 
  geom_boxplot() +
  geom_quasirandom()+
  geom_signif(comparisons = list(c("0", "1"), c("0", "2"), c("0", "3+"), 
                                 c("1", "2"), c("1", "3+"), c("2", "3+")), 
              y_position = c(9,9.5,10,8, 8.5, 9), test = t.test, 
              map_signif_level=T)+
  facet_wrap(MetricL~., scales = "free_y")+ theme_bw()+
 # geom_text(data = pairs, aes(y = Y, label = Group))+
  xlab("Drought Length (Years)")+ ylab("log FMWT index")+ ylim(3,11)+
  scale_fill_brewer(palette = "Dark2", name = NULL, guide = NULL)
sh
```


```{r Fig3b, warn = FALSE }
# Srcrr = ggplot(SR,  aes(x = as.factor(DroughtYear), y = Value, fill = as.factor(DroughtYear))) + 
#   geom_boxplot() +
#   geom_quasirandom()+
#   geom_signif(comparisons = list(c("0", "1"), c("0", "2"), c("0", "3+"), 
#                                  c("1", "2"), c("1", "3+"), c("2", "3+")), 
#               y_position = c(4,4.5,5,4.2, 3.7), test = t.test, 
#               map_signif_level=T)+
#   facet_wrap(MetricL~., scales = "free_y")+ theme_bw()+
#  # geom_text(data = pairs, aes(y = Y, label = Group))+
#   xlab("Return Year")+ ylab("CRR")+ ylim(0,5.5)+
#   scale_fill_brewer(palette = "Dark2", name = NULL, guide = NULL)
# 
#   Srcrr
  
    SR = mutate(SR, MetricL = "Spring Run Chinook") %>%
    filter(!is.na(DroughtYearSalmon))
  Srcrr2 = ggplot(SR,  aes(x = as.factor(DroughtYearSalmon), y = Value, fill = as.factor(DroughtYearSalmon))) + 
  geom_boxplot() +
  geom_quasirandom()+
  geom_signif(comparisons = list(c("0", "1"), c("0", "2"), c("0", "3+"), 
                                 c("1", "2"), c("1", "3+"), c("2", "3+")), 
              y_position = c(4,4.5,5,4.2, 3.7, 3.4), test = t.test, 
              map_signif_level=T)+
  facet_wrap(MetricL~., scales = "free_y")+ theme_bw()+
 # geom_text(data = pairs, aes(y = Y, label = Group))+
  xlab("Drought Length (Migration Year)")+ ylab("CRR")+ ylim(0,5.5)+
  scale_fill_brewer(palette = "Dark2", name = NULL, guide = NULL)

  Srcrr2
```


```{r Fig3c, warn = FALSE }
# Wrcrr = ggplot(WR,  aes(x = as.factor(DroughtYear), y = Value, fill = as.factor(DroughtYear))) + 
#   geom_boxplot() +
#   geom_quasirandom()+
#       geom_signif(comparisons = list(c("0", "1"), c("0", "2"), c("0", "3+"), 
#                                  c("1", "2"), c("1", "3+"), c("2", "3+")), 
#               y_position = c(4,4.5,5,4.2, 3.7), test = t.test, 
#               map_signif_level=T)+
#   facet_wrap(MetricL~., scales = "free_y")+ theme_bw()+
#  # geom_text(data = pairs, aes(y = Y, label = Group))+
#   xlab("Return Year")+ ylab("CRR")+ #ylim(0,5.5)+
#   scale_fill_brewer(palette = "Dark2", name = NULL, guide = NULL)
# 
#   Wrcrr
  
  WR = mutate(WR, MetricL = "Winter Run Chinook") %>%
    filter(!is.na(DroughtYearSalmon))
  Wrcrr2 = ggplot(WR,  aes(x = as.factor(DroughtYearSalmon), y = Value, fill = as.factor(DroughtYearSalmon))) + 
  geom_boxplot() +
  geom_quasirandom()+
    geom_signif(comparisons = list(c("0", "1"), c("0", "2"), c("0", "3+"), 
                                 c("1", "2"), c("1", "3+"), c("2", "3+")), 
              y_position = c(4,4.5,5,4.2, 3.7, 3.4), test = t.test, 
              map_signif_level=T)+
  facet_wrap(MetricL~., scales = "free_y")+ theme_bw()+
 # geom_text(data = pairs, aes(y = Y, label = Group))+
  xlab("Drought Length (Migration Year)")+ ylab("CRR")+ #ylim(0,5.5)+
  scale_fill_brewer(palette = "Dark2", name = NULL, guide = NULL)

  
  Wrcrr2
```


```{r Fig3d, warn = FALSE }
# 
#  Frcrr = ggplot(FR,  aes(x = as.factor(DroughtYear), y = Value, fill = as.factor(DroughtYear))) + 
#   geom_boxplot() +
#   geom_quasirandom()+
#   geom_signif(comparisons = list(c("0", "1"), c("0", "2"), c("0", "3+"), 
#                                  c("1", "2"), c("1", "3+"), c("2", "3+")), 
#               y_position = c(4,4.5,5,4.2, 3.7), test = t.test, 
#               map_signif_level=T)+
#   facet_wrap(MetricL~., scales = "free_y")+ theme_bw()+
#  # geom_text(data = pairs, aes(y = Y, label = Group))+
#   xlab("return year")+ ylab(NULL)+ ylim(0,5.5)+
#   scale_fill_brewer(palette = "Dark2", name = NULL, guide = NULL)
# 
#   Frcrr
  FR = mutate(WR, MetricL = "Fall Run Chinook") %>%
    filter(!is.na(DroughtYearSalmon))
 Frcrr2 = ggplot(FR,  aes(x = as.factor(DroughtYearSalmon), y = Value, fill = as.factor(DroughtYearSalmon))) + 
  geom_boxplot() +
  geom_quasirandom()+
  geom_signif(comparisons = list(c("0", "1"), c("0", "2"), c("0", "3+"), 
                                 c("1", "2"), c("1", "3+"), c("2", "3+")), 
              y_position = c(4,4.5,5,4.2, 3.7, 3.4), test = t.test, 
              map_signif_level=T)+
  facet_wrap(MetricL~., scales = "free_y")+ theme_bw()+
 # geom_text(data = pairs, aes(y = Y, label = Group))+
  xlab("Drought Length (Migration Year)")+ ylab("CRR")+ ylim(0,5.5)+
  scale_fill_brewer(palette = "Dark2", name = NULL, guide = NULL)

  Frcrr2
```


```{r Fig3d, warn = FALSE }
ds = ggplot(filter(FishDrought, Metric == "logDS"),  aes(x = as.factor(DroughtYear), y = Value, fill = as.factor(DroughtYear))) + 
  geom_boxplot() +
  geom_quasirandom()+
  geom_signif(comparisons = list(c("0", "3+")), annotations = "All comparisons NS")+
  facet_wrap(MetricL~., scales = "free_y")+ theme_bw()+
 # geom_text(data = pairs, aes(y = Y, label = Group))+
  xlab("Drought Length (Years)")+ ylab("log+1 FMWT index")+ ylim(0,9)+
  scale_fill_brewer(palette = "Dark2", name = NULL, guide = NULL)+

  theme(legend.position = NULL)

  ds
```


```{r Fig3d, warn = FALSE }
sb = ggplot(filter(FishDrought, Metric == "logSB", !is.na(Value)),  aes(x = as.factor(DroughtYear), y = Value, fill = as.factor(DroughtYear))) + 
  geom_boxplot() +
  geom_quasirandom()+
  geom_signif(comparisons = list(c("0", "1"), c("0", "2"), c("0", "3+"), 
                                 c("1", "2"), c("1", "3+"), c("2", "3+")), 
              y_position = c(10,11,12,8.5, 9.5, 8), test = t.test, 
              map_signif_level=T)+
  facet_wrap(MetricL~., scales = "free_y")+ theme_bw()+
 # geom_text(data = pairs, aes(y = Y, label = Group))+
  xlab("Drought Length (Years)")+ ylab("log FMWT index")+ ylim(0,12.5)+
  scale_fill_brewer(palette = "Dark2", name = NULL, labels = c("0, Wet years", "1, First/only dry year",
                                                                            "2, Second dry year", "3+ Third or more dry year"))+
    theme(legend.position = "bottom")

  sb
  sb2 = sb +  scale_fill_brewer(palette = "Dark2", guide = NULL)
```


```{r Fig3d, warn = FALSE }
TFS =  ggplot(filter(FishDrought, Metric == "logTFS"),  aes(x = as.factor(DroughtYear), y = Value, fill = as.factor(DroughtYear))) + 
  geom_boxplot() +
  geom_quasirandom()+
  geom_signif(comparisons = list(c("0", "3+")), annotations = "All comparisons NS")+
  facet_wrap(MetricL~., scales = "free_y")+ theme_bw()+
 # geom_text(data = pairs, aes(y = Y, label = Group))+
  xlab("Drought Length (Years)")+ ylab("log FMWT index")+ ylim(0,12)+
  scale_fill_brewer(palette = "Dark2", name = NULL, guide = NULL)+

  theme(legend.position = NULL)

  TFS
```


```{r Fig3d, warn = FALSE }
#print just the legend
legend <- cowplot::get_legend(sb)

plots = grid.arrange(ds, LF, sh,TFS, sb2, Wrcrr2,Frcrr2, Srcrr2, legend,  
                     heights = c(2,2,2,2,.5),
                     layout_matrix = rbind(c(1, 2), c(3,4), c(5,6), c(7,8), c(9,9)))

plotsnosalmon = grid.arrange(ds, LF, sh,TFS, sb2, legend,  
                     nrow =4, ncol = 2,heights = c(2,2,2,2))
  
####THIS IS FIGURE 3
ggsave(plot = plots,"plots/FishDrought.tiff", device = "tiff", width = 8, height = 12, units = "in")


#now do just the salmon

salmon = grid.arrange(Wrcrr2,Frcrr2,  Srcrr2,legend,
                     nrow =2, ncol = 2,heights = c(2,2))

  
ggsave(plot = salmon,"plots/SalmonDrought.tiff", device = "tiff", width = 8, height = 12, units = "in")


```
Now i need a table with all the statistical results

```{r}

library(broom)

FD1t = tidy(FD1) %>%
  mutate(Species = "Striped Bass")

LFS1t = tidy(LFS1) %>%
  mutate(Species = "Longfin Smelt")

shad1t = tidy(shad1) %>%
  mutate(Species = "American Shad")


Deltas = tidy(DS1) %>%
  mutate(Species = "Delta Smelt")

Spring = tidy(SRcrr) %>%
  mutate(Species = "Spring Run CRR")
Winter = tidy(WRcrr) %>%
  mutate(Species = "Winter Run CRR")

Fall = tidy(FRcrr) %>%
  mutate(Species = "Fall Run CRR")

Threadfin = tidy(TFshad1) %>%
  mutate(Species = "Threadfin Shad")

models = bind_rows(FD1t, LFS1t, shad1t, Deltas, Threadfin, Spring, Winter, Fall)
write.csv(models, "outputs/FishMultiYear.csv", row.names = F)
```
Plots for white papers

```{r, warn = FALSE}
library(ggbeeswarm)
FishDrought2 = filter(FishDrought, Metric != "Salmon CRR") %>%
  mutate(Whitepaper = factor(Whitepaper, levels = c("Critical", "Dry", "Below Normal", "Above Normal", "Wet", "2020", "2021", "2022"),
                             labels = c("Critical", "Dry", "Below\nNormal", "Above\nNormal", "Wet", "2020", "2021", "2022")))
  

ggplot(FishDrought2,  aes(x = Whitepaper, y = Value, fill = Yr_type)) + 
  geom_boxplot() +
  drt_color_pal_yrtype()+
  geom_quasirandom()+
  facet_wrap(MetricL~., scales = "free_y", nrow = 5)+ theme_bw()+
   ylab("log FMWT Index")+xlab(NULL)+
  theme(legend.position = "none")

ggsave("plots/whitepaper/fish.tiff", device = "tiff", height = 8, width = 5)

```