results_summary.qmd

---
title: "Synchrony Paper: Results"
author: "R. Santos and W.R. James"
date: "1/30/2024"
format: 
  html:
    toc: true
    theme: yeti
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(collapse = T, cache = T)
```

## Content: 
This document includes the analyses for

1.  Q1: How the inter-individual variability in space use varies across seasons and changed over the years?
      + Eadj ~ s(Seasons) + s(Year)
2.  Q2: How the inter-individual variability in space use respond to hydrological parameters
      + Eadj ~ s(water stage) + s(days below 30 cm)
3.  Q3: Does the inter-individual variability in space use translate into trophic niche changes?   


R script: [github](https://github.com/CoastalFishScience/spatsim/blob/main/script/stats_analyses.R)
      

##  Packages
Packages with functions used for data processing and analyses. 
```{r}
library(tidyverse)
library(ggpubr)

#Models libraries
library(mgcv)
library(glmmTMB)

#Model check libraries
library(visreg)
library(DHARMa)
library(performance)
library(lsmeans)
library(emmeans)
library(MuMIn)
library(ggeffects)
```


## Data processing

Data with Eadj metric, date, hydrology based seasons and year, and hydrological metrics. Metrics are based on monthly observations and divided by water year (May-April). 

For script used to estimate Eadj values please go to:
R script: [Eadj](https://github.com/CoastalFishScience/spatsim/blob/main/script/SpatSimEadj_2024_MW.R)

For script used to estimate trophic niche size got to:
R script: [MixingModels](https://github.com/CoastalFishScience/spatsim/blob/main/script/SpatSimMixingModels_2024_MW.R)

R script: [Hypervolume](https://github.com/CoastalFishScience/spatsim/blob/main/script/SpatSimHyperVolumes_2024_MW.R)

```{r}
#all data
sim = read_csv("./data/spat_sim_allthegoods_01_15_2023.csv")
glimpse(sim)

#sample size info
nsize = read_csv("./data/snook_sample_size_eadj.csv") |>  
      dplyr::rename(Year.Month = fYear.Month, eadj_n = n)
glimpse(nsize)

sim_all <- left_join(sim, nsize, by = "Year.Month")
glimpse(sim_all)
```


## Q1 analyses, results and plots 
Q1: How the inter-individual variability in space use varies across seasons and changed over the years?

We used GAMs to understand and characterize the seasonal cycle in inter-individual variability in space use and the overall trend across the years

### Model Selection:
```{r}

#Both smoothers for water months and year
q1.m1 <- gam(Eadj ~ s(wMonth, bs = "cc", k = 10) + s(wYear), 
          data = sim_all,
          family = betar(link = "logit"), 
          method = "REML")

#Smoother only  for water months
q1.m2 <- gam(Eadj ~ s(wMonth, bs = "cc", k = 10) + wYear, 
             data = sim_all,
             family = betar(link = "logit"), 
             method = "REML")

#GLM (no smoothers)
q1.m3 <- glmmTMB(Eadj ~ wMonth + wYear, 
             data = sim_all,
             family = beta_family(link = "logit"))

#Compare performance between models
compare_performance(q1.m1, q1.m2, q1.m3)
```

Based on AICc the model with both smoothers (q1.m1) was selected as the best model structure

### Model (q1.m1) summary
Summary of the coefficients estimates and goodness of fit


```{r}
summary(q1.m1) 
```

### Plots of fitted model  
```{r}
#Model output - seasonal trend
q1.m1.vis = visreg(q1.m1, type = "conditional", scale = "response")
#extracting fitted model from visreg to allow plotting in ggplot
q1.m1.month<-q1.m1.vis[[1]]$fit

seasonal_trend<-ggplot(q1.m1.month, aes(wMonth, visregFit))+ 
      geom_ribbon(aes(ymin = visregLwr, ymax = visregUpr), color = "grey60", alpha = .2)+
      geom_line(linewidth = 2, colour= "black") + theme_bw()+
      geom_line(linetype = 2, colour = "black", aes(y = visregLwr))+
      geom_line(linetype = 2, colour = "black", aes(y = visregUpr)) +
      labs(x = "Water Year Month (May-April)", y = "Eadj", title = "Seasonal trend")+
      scale_x_continuous(breaks = c(2, 4, 6, 8, 10, 12)) +
      ylim(0.3,0.8) +
      theme(axis.text = element_text(size = 14, face = "bold", colour = "black"), 
            axis.title = element_text(size = 15, face = "bold", colour = "black"), 
            plot.title = element_text(size = 16, face = "bold", colour = "black", hjust = 0.5), 
            panel.grid.major = element_blank(),
            axis.line = element_line(colour = "black"),
            panel.grid.minor = element_blank(),
            panel.border = element_blank(),
            panel.background = element_blank()) 

#Model ouptut - Yearly variability 
q1.m1.year<-q1.m1.vis[[2]]$fit

yearly_vars<-ggplot(q1.m1.year, aes(wYear, visregFit))+ 
      geom_ribbon(aes(ymin = visregLwr, ymax = visregUpr), color = "grey60", alpha = .2)+
      geom_line(size = 2, colour= "black") + theme_bw()+
      geom_line(linetype = 2, colour = "black", aes(y = visregLwr))+
      geom_line(linetype = 2, colour = "black", aes(y = visregUpr)) +
      labs(x = "Water Year (May-April)", y = "Eadj", title = "Inter-annual trend")+
      scale_x_continuous(breaks = c(2012, 2014, 2016, 2018, 2020, 2022)) +
      ylim(0.3,0.8)+
      theme(axis.text = element_text(size = 14, face = "bold", colour = "black"), 
            axis.title = element_text(size = 15, face = "bold", colour = "black"), 
            plot.title = element_text(size = 16, face = "bold", colour = "black", hjust = 0.5), 
            panel.grid.major = element_blank(),
            axis.line = element_line(colour = "black"),
            panel.grid.minor = element_blank(),
            panel.border = element_blank(),
            panel.background = element_blank())


#Plot of fitted model q1.m1
ggarrange(seasonal_trend, yearly_vars,
          labels = c('a)','b)'),
          ncol = 1, vjust = 1, align = "v")
```


## Q2 analyses, results and plots 

Q2: How the inter-individual variability in space use respond to hydrological parameters 

We used GLMs to understand and characterize the response in inter-individual variability in space use to hydrological conditions

### Model Selection: Due to collinearity only single variable models were performed
```{r}

#Renaming variables for model
sim_all = sim_all |> 
      rename(mean_stage = monthly_mean_stage_cm,
             mean_Lstage = monthly_mean_stage_cm_prev)

#GLM and GAM using mean_stage
q2.glm.stage = glmmTMB(Eadj ~ mean_stage, 
                    data = sim_all,
                    family = beta_family(link = "logit"))

q2.gam.stage = gam(Eadj ~ s(mean_stage), 
                   data = sim_all,
                   family = betar(link = "logit"), 
                   method = "REML")



#GLM and GAM using days below 30 cm
q2.glm.days = glmmTMB(Eadj ~ daysbelow30, 
                        data = sim_all,
                        family = beta_family(link = "logit"))

q2.gam.days = gam(Eadj ~ s(daysbelow30), 
                  data = sim_all,
                  family = betar(link = "logit"), 
                  method = "REML")




#GLM and GAM using lag of mean stage
q2.gam.Lstage = gam(Eadj ~ s(mean_Lstage), 
                   data = sim_all,
                   family = betar(link = "logit"), 
                   method = "REML")

q2.glm.Lstage = glmmTMB(Eadj ~ mean_Lstage, 
                       data = sim_all,
                       family = beta_family(link = "logit"))

#Best = gam
compare_performance(q2.glm.stage, q2.gam.stage)
compare_performance(q2.glm.Lstage, q2.gam.Lstage)  
compare_performance(q2.glm.days, q2.gam.days) 
```

Based on AICc GAMs provided the best model structure. However, also based on AICc, we could not determined a single best model - i.e., single variable models were indistinguishable based on AICc. Based on these results we decided to report the results both of the single variable models

### Model (q2.gam.stage, q2.gam.days) summary
Summary of the coefficients estimates and goodness of fit


```{r}
summary(q2.gam.stage)
summary(q2.gam.days)
```

### Plots of fitted model  
```{r}
#daysbelow30 effects
q2.days.vis = visreg(q2.gam.days, type = "conditional", scale = "response")
q2.days.fit<-q2.days.vis$fit

days_vs_E<-ggplot(q2.days.fit, aes(daysbelow30, visregFit))+
      geom_ribbon(aes(ymin = visregLwr, ymax = visregUpr), color = "grey60", alpha = .2)+
      geom_line(size = 2, colour= "black") + theme_bw()+
      geom_line(linetype = 2, colour = "black", aes(y = visregLwr))+
      geom_line(linetype = 2, colour = "black", aes(y = visregUpr)) +
      labs(x = "Days with stage below 30 cm", y = "Eadj")+
      scale_x_continuous(breaks = c(0, 5, 10, 15, 20, 25, 30)) +
      ylim(0.3,0.8)+
      theme(axis.text = element_text(size = 14, face = "bold", colour = "black"),
            axis.title = element_text(size = 16, face = "bold", colour = "black"),
            plot.title = element_text(size = 16, face = "bold", colour = "black"),
            panel.grid.major = element_blank(),
            axis.line = element_line(colour = "black"),
            panel.grid.minor = element_blank(),
            panel.border = element_blank(),
            panel.background = element_blank())
      
      
#water stage effects
q2.stage.vis = visreg(q2.gam.stage, type = "conditional", scale = "response")
q2.stage.fit<-q2.stage.vis$fit

stage_vs_E<-ggplot(q2.stage.fit, aes(mean_stage, visregFit))+
      geom_ribbon(aes(ymin = visregLwr, ymax = visregUpr), color = "grey60", alpha = .2)+
      geom_line(size = 2, colour= "black") + theme_bw()+
      geom_line(linetype = 2, colour = "black", aes(y = visregLwr))+
      geom_line(linetype = 2, colour = "black", aes(y = visregUpr)) +
      labs(x = "River Stage (cm)", y = "Eadj")+
      scale_x_continuous(breaks = c(0, 5, 10, 15, 20, 25, 30)) +
      ylim(0.3,0.8)+
      theme(axis.text = element_text(size = 14, face = "bold", colour = "black"),
            axis.title = element_text(size = 16, face = "bold", colour = "black"),
            plot.title = element_text(size = 16, face = "bold", colour = "black"),
            panel.grid.major = element_blank(),
            axis.line = element_line(colour = "black"),
            panel.grid.minor = element_blank(),
            panel.border = element_blank(),
            panel.background = element_blank())

#Plot of fitted model q1.m1
ggarrange(stage_vs_E, days_vs_E,
          labels = c('a)','b)'),
          ncol = 1, vjust = 1, align = "v")

```

## Q3 analyses, results and plots 
Q3: Does the inter-individual variability in space use translate into trophic niche changes?   

We used pearson correlation to assess the assoication of inter-individual variability in space use with trophic niche size

### Preparing dataset

```{r}
#data wraggling to correlate dry montly Eadj and HV
dat_summary <- sim |> 
      group_by(wYear, Season) |> 
      summarise(mean_e = mean(Eadj, na.rm = "TRUE"),
                min_e = min(Eadj, na.rm = "TRUE"),
                max_e = max(Eadj, na.rm = "TRUE"),
                hv_size = mean(hv_size, na.rm = "TRUE"),
                hv_n = mean(hv_n)) |> 
      filter(Season == "Dry")
```

### Correlation analysis
```{r}

#All data
cor.test(dat_summary$hv_size, dat_summary$mean_e, method = "pearson")

#Without large outlier
dat_summary2 = filter(dat_summary, hv_size < 200) 
cor.test(dat_summary2$hv_size, dat_summary2$mean_e, method = "pearson")
```

Based on correlation analysis, only significant association without the outlier


### Plots of fitted model  
```{r}
corr_all = ggplot(dat_summary, aes(x = mean_e, y = hv_size)) + 
      geom_point() +
      geom_smooth(method = "lm", se = TRUE, color = "black") +
      labs(x = "Eadj Dry Season Mean", 
           y = "Niche Volume") +
      stat_cor(p.accuracy = 0.01, r.accuracy = 0.01, label.x.npc = .5)+
      xlim(0.4,0.7) +
      # ylim(0, max(dat_summary$hv_size))+
      theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
            panel.background = element_blank(), axis.line = element_line(colour = "black")) + 
      theme(plot.title = element_text(hjust = 0.5)) +
      theme(plot.title = element_text(size=14, face="bold", color = "black")) +
      theme(axis.text = element_text(size=16,face="bold", color = "black")) +
      theme(axis.text.x = element_text(size=16,face="bold", color = "black")) +
      theme(axis.text.y = element_text(size=16,face="bold", color = "black")) +
      theme(axis.title = element_text(size=16,face="bold", color = "black")) +
      theme(legend.title = element_blank()) +
      theme(legend.text = element_text(size=16, face="bold", color = "black")) +
      theme(legend.position = c(0.9, 0.9))

corr_NOoutlier = ggplot(filter(dat_summary, hv_size < 200), aes(x = mean_e, y = hv_size)) + 
      geom_point() +
      geom_smooth(method = "lm", se = TRUE, color = "black") +
      labs(x = "Eadj Dry Season Mean", 
           y = "Niche Volume") +
      stat_cor(p.accuracy = 0.01, r.accuracy = 0.01, label.x.npc = .5)+
      xlim(0.4,0.7) +
      # ylim()
      theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
            panel.background = element_blank(), axis.line = element_line(colour = "black")) + 
      theme(plot.title = element_text(hjust = 0.5)) +
      theme(plot.title = element_text(size=14, face="bold", color = "black")) +
      theme(axis.text = element_text(size=16,face="bold", color = "black")) +
      theme(axis.text.x = element_text(size=16,face="bold", color = "black")) +
      theme(axis.text.y = element_text(size=16,face="bold", color = "black")) +
      theme(axis.title = element_text(size=16,face="bold", color = "black")) +
      theme(legend.title = element_blank()) +
      theme(legend.text = element_text(size=16, face="bold", color = "black")) +
      theme(legend.position = c(0.9, 0.9))

#Plot of fitted model q1.m1
ggarrange(corr_all, corr_NOoutlier,
          labels = c('a)','b)'),
          ncol = 1, vjust = 1, align = "v")

```