---
title: "Kastl_2022_modeling"
author: "Brian Kastl"
date: "1/29/2021"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


```{r}
library(data.table)
library(readxl)
library(readr)
library(devtools)
library(rgdal)
library(foreign)
library(tidyverse)
library(dplyr)
library(sp)
library(sf)
library(ggplot2)
library(stringr)
library(tidyr)
library(stringr)
library(lubridate)
library(ggpmisc)
library(spatstat)
library(GGally)
library(ggpmisc)
```


```{r}
library("Hmisc")

# Read in By Stream Year data from Vis_Stats53.Rmd
Output_StYr_all <- read_csv('../output/Kastl_2022_modeling_input.csv')
```


FILTER OUT STREAM-YEARS WITH < 3 YEARS OF DATA (for Lmer models)
```{r}
Output_StYr_all <- Output_StYr_all %>%
  filter(!(Stream == 'Green V. LST'))
```


```{r}
mean_start <- mean(Output_StYr_all$Five_pct)
mean_start

#St dev: divide by number of secs per day
sd_start <- sd(Output_StYr_all$Five_pct)/86400
sd_start

mean_med <- mean(Output_StYr_all$Fifty_pct)
mean_med

#St dev: divide by number of secs per day
sd_med <- sd(Output_StYr_all$Fifty_pct)/86400
sd_med

mean_end <- mean(Output_StYr_all$NinetyFive_pct)
mean_end

#St dev: divide by number of secs per day
sd_end <- sd(Output_StYr_all$NinetyFive_pct)/86400
sd_end

mean_dur <- mean(Output_StYr_all$Dur_5_95)
mean_dur

sd_dur <- sd(Output_StYr_all$Dur_5_95)
sd_dur
```


```{r}
# Select only needed columns
Output_StYr_all_trim <- Output_StYr_all %>%
  select(Degree_Days_sum_MarApr, MarJun_Trib_Runoff_mmPerDay, MarApr_Trib_Runoff_mmPerDay) %>%
  # delete all rows with any NA in any column (to make it the 42 stream-years, used in modeling)
  drop_na()

# Correlation coefficients (using only rows without any NAs) (for some reason, these values are slightly diff than further below [maybe b/c "complete.obs"])
pearson_all <- cor(Output_StYr_all_trim, method = "pearson", use = "complete.obs") %>% round(2)
spearman_all <- cor(Output_StYr_all_trim, method = "spearman", use = "complete.obs") %>% round(2)

# Correlation matrix with significance levels (p-value)
pearson_all2 <- rcorr(as.matrix(Output_StYr_all_trim), type = "pearson")
spearman_all2 <- rcorr(as.matrix(Output_StYr_all_trim), type = "spearman")


# flattenCorrMatrix
# cormat : matrix of the correlation coefficients
# pmat : matrix of the correlation p-values
flattenCorrMatrix <- function(cormat, pmat) {
  ut <- upper.tri(cormat)
  data.frame(
    row = rownames(cormat)[row(cormat)[ut]],
    column = rownames(cormat)[col(cormat)[ut]],
    cor  =(cormat)[ut],
    p = pmat[ut]
    )
}

Flat_pearson_all <- flattenCorrMatrix(pearson_all2$r, pearson_all2$P)
Flat_pearson_all

Flat_spearman_all <- flattenCorrMatrix(spearman_all2$r, spearman_all2$P)
Flat_spearman_all
```


Selecting the needed variables for the LMER
```{r}
Output_StYr_all_lmer <- Output_StYr_all %>%
    select(Stream, Year, Dur_5_95, Five_pct, NinetyFive_pct, Degree_Days_sum_MarApr, MarJun_Trib_Runoff_mmPerDay, MarApr_Trib_Runoff_mmPerDay) %>%
  # To address warning: "Some predictor variables are on very different scales: consider rescaling"
  mutate(Degree_Days_sum_scaled_MarApr = Degree_Days_sum_MarApr / 100)
```
  

Convert date to numeric for 5pct and 95pct.
In two chunks below, running the lmer I had non-convergence error, so I divide by 1e9
```{r}
# as.numeric needed to make model run, which converts date to "seconds after a date in history"
# divide by the number of seconds per day to make the units Five_pct and NinetyFive_pct into "days"
Output_StYr_all_lmer$Five_pct <- as.numeric(Output_StYr_all_lmer$Five_pct)/86400

Output_StYr_all_lmer$NinetyFive_pct <- as.numeric(Output_StYr_all_lmer$NinetyFive_pct)/86400
```

 
Since I can have different sample sizes for models with different response variables (AIC values are only compared among models with the same response variable), create data for which NA rows are only dropped for explanatory variables that are used for the response variable
```{r} 
Output_StYr_all_lmer_MarApr <- Output_StYr_all_lmer %>%
  select(Stream, Year, Dur_5_95, Five_pct, NinetyFive_pct, Degree_Days_sum_MarApr, Degree_Days_sum_scaled_MarApr, MarApr_Trib_Runoff_mmPerDay, MarJun_Trib_Runoff_mmPerDay) %>%
# Eliminate any row that has an NA value in any column
  drop_na()
```


WRITE OUT TABLES
```{r}
# WRITE OUT SUPPLEMENTAL TABLE THAT HAS THE FINAL DATA WE USE TO RUN THE MODEL
Output_StYr_all_lmer_select <- Output_StYr_all_lmer_MarApr %>%
  select(Stream, Year, Dur_5_95, Five_pct, NinetyFive_pct, Degree_Days_sum_MarApr, Degree_Days_sum_scaled_MarApr, MarApr_Trib_Runoff_mmPerDay, MarJun_Trib_Runoff_mmPerDay)

Output_StYr_all_lmer_select$Five_pct <- format(as.Date(Output_StYr_all_lmer_select$Five_pct, origin = "1970-01-01 00:00:00"), "%m/%d")
Output_StYr_all_lmer_select$NinetyFive_pct <- format(as.Date(Output_StYr_all_lmer_select$NinetyFive_pct, origin = "1970-01-01 00:00:00"), "%m/%d")

write.csv(Output_StYr_all_lmer_select, '../output/Tables_appendix/Kastl_2022_model_input.csv')


# WRITE OUT MAIN PAPER TABLE OF STUDY AREA CHARACTERISTICS
# Using earlier, pre drop.na() version because Ted wants to show the number of stream-years of outmigration data pre-filtering
Output_StYr_all_byStrm <- Output_StYr_all %>%
  group_by(Stream) %>%
  summarise(outmigration_years = length(which(!is.na(Five_pct))),
            degree_days_MarApr_years = length(which(!is.na(Degree_Days_sum_MarApr))),
            degree_days_MarApr_min = min(Degree_Days_sum_MarApr, na.rm=TRUE),
            degree_days_MarApr_max = max(Degree_Days_sum_MarApr, na.rm=TRUE),
            runoff_mmPerDay_MarApr_min = min(MarApr_Trib_Runoff_mmPerDay),
            runoff_mmPerDay_MarApr_max = max(MarApr_Trib_Runoff_mmPerDay),
            runoff_mmPerDay_MarJun_min = min(MarJun_Trib_Runoff_mmPerDay),
            runoff_mmPerDay_MarJun_max = max(MarJun_Trib_Runoff_mmPerDay))

DA_Precip2 <- read.csv('../data/Kastl_2022_DrainageArea_Precip.csv')

Study_area_table <- left_join(Output_StYr_all_byStrm, DA_Precip2, by = "Stream") %>%
  select(Stream, outmigration_years, degree_days_MarApr_years, degree_days_MarApr_min, degree_days_MarApr_max, runoff_mmPerDay_MarApr_min, runoff_mmPerDay_MarApr_max, runoff_mmPerDay_MarJun_min, runoff_mmPerDay_MarJun_max, Arc_Drainage_sqkm, Precip_mm_mean2011_2020)

write.csv(Study_area_table, '../output/Tables_appendix/Kastl_study_area_table.csv')
```


LMER: Response variable = Five_pct (USE DATA: Output_StYr_all_lmer_MarApr)
Beginner's guide: https://rstudio-pubs-static.s3.amazonaws.com/63556_e35cc7e2dfb54a5bb551f3fa4b3ec4ae.html
```{r}
library(lme4)
library(lmerTest)
library(psycho)
library(robustlmm)
require(stats)

# As interactive ("*") with MarApr Runoff
L5_MarAprRunoff_MarAprDegree_Int <- lmer(Five_pct ~ MarApr_Trib_Runoff_mmPerDay * Degree_Days_sum_scaled_MarApr + (1|Stream), data = Output_StYr_all_lmer_MarApr, na.action = na.fail, REML = FALSE)

# As additive ("+") with MarApr Runoff
L5_MarAprRunoff_MarAprDegree_Ad <- lmer(Five_pct ~ MarApr_Trib_Runoff_mmPerDay + Degree_Days_sum_scaled_MarApr + (1|Stream), data = Output_StYr_all_lmer_MarApr, na.action = na.fail, REML = FALSE)

# Date = Output_StYr_all_lmer b/c able to use 47 stream-years
L5_MarAprRunoff_only <- lmer(Five_pct ~ MarApr_Trib_Runoff_mmPerDay + (1|Stream), data = Output_StYr_all_lmer_MarApr, na.action = na.fail, REML = FALSE)

L5_MarAprDegree_only <- lmer(Five_pct ~ Degree_Days_sum_scaled_MarApr + (1|Stream), data = Output_StYr_all_lmer_MarApr, na.action = na.fail, REML = FALSE)
```

LMER: Response variable = NinetyFive_pct (USE DATA: Output_StYr_all_lmer_MarJun)
```{r}
# As interactive ("*") with MarJun Runoff
L95_MarJunRunoff_MarAprDegree_Int <- lmer(NinetyFive_pct ~ MarJun_Trib_Runoff_mmPerDay * Degree_Days_sum_scaled_MarApr + (1|Stream), data = Output_StYr_all_lmer_MarApr, na.action = na.fail, REML = FALSE)

# As additive ("+") with MarJun Runoff + MarJun degree days
L95_MarJunRunoff_MarAprDegree_Ad <- lmer(NinetyFive_pct ~ MarJun_Trib_Runoff_mmPerDay + Degree_Days_sum_scaled_MarApr + (1|Stream), data = Output_StYr_all_lmer_MarApr, na.action = na.fail, REML = FALSE)

# Date = Output_StYr_all_lmer b/c able to use 47 stream-years
L95_MarJunRunoff_only <- lmer(NinetyFive_pct ~ MarJun_Trib_Runoff_mmPerDay + (1|Stream), data = Output_StYr_all_lmer_MarApr, na.action = na.fail, REML = FALSE)

L95_MarAprDegree_only <- lmer(NinetyFive_pct ~ Degree_Days_sum_scaled_MarApr + (1|Stream), data = Output_StYr_all_lmer_MarApr, na.action = na.fail, REML = FALSE)
```


LMER: Response variable = Dur_5_95  (USE DATA: Output_StYr_all_lmer_MarJun)
```{r}
L_MarJunRunoff_MarAprDegree_Int <- lmer(Dur_5_95 ~ MarJun_Trib_Runoff_mmPerDay * Degree_Days_sum_scaled_MarApr + (1|Stream), data = Output_StYr_all_lmer_MarApr, na.action = na.fail, REML = FALSE)

L_MarJunRunoff_MarAprDegree_Ad <- lmer(Dur_5_95 ~ MarJun_Trib_Runoff_mmPerDay + Degree_Days_sum_scaled_MarApr + (1|Stream), data = Output_StYr_all_lmer_MarApr, na.action = na.fail, REML = FALSE)

# As only MarJun RUNOFF for explanatory variable
# USING Output_StYr_all_lmer (47 stream-years), INSTEAD OF Output_StYr_all_lmer_MarApr (42 stream-years), BECAUSE WE HAVE RUNOFF ONLY AND CAN USE 5 EXTRA DAYS THIS WAY (TED SUGGESTED)
L_MarJunRunoff_only <- lmer(Dur_5_95 ~ MarJun_Trib_Runoff_mmPerDay + (1|Stream), data = Output_StYr_all_lmer_MarApr, na.action = na.fail, REML = FALSE)

L_MarAprDegree_only <- lmer(Dur_5_95 ~ Degree_Days_sum_scaled_MarApr + (1|Stream), data = Output_StYr_all_lmer_MarApr, na.action = na.fail, REML = FALSE)
```


COMPARE STRENGTH OF MODELS
AIC Model comparison (the lowest score is the strongest)
```{r}
# Five_pct
Gabe_anova_5 <- anova(L5_MarAprRunoff_MarAprDegree_Int, L5_MarAprRunoff_MarAprDegree_Ad, L5_MarAprDegree_only, L5_MarAprRunoff_only) %>% 
  rownames_to_column() %>%
  rename("Model" = rowname)

# NinetyFive_pct
Gabe_anova_95 <- anova(L95_MarJunRunoff_MarAprDegree_Int, L95_MarJunRunoff_MarAprDegree_Ad, L95_MarJunRunoff_only, L95_MarAprDegree_only) %>%
  rownames_to_column() %>%
  rename("Model" = rowname)

# 5_95
Gabe_anova <- anova(L_MarJunRunoff_MarAprDegree_Ad, L_MarJunRunoff_only, L_MarAprDegree_only) %>%
  rownames_to_column() %>%
  rename("Model" = rowname)
```


```{r}
AIC_table <- rbind(Gabe_anova_5, Gabe_anova_95, Gabe_anova)
```


5_pct: R squared values for fixed effects you can use the r.squaredGLMM() command from the ‘MuMIn package. It shows the marginal and the conditional R².The marginal R2 is for only the fixed effects and conditional R2 shows the variance for the entire model (including fixed and random effects). 
I manually typed in the r.squaredGLMM into Summary_p_values_Rsq.csv.
```{r}
library(MuMIn)
library(gdata)

# For top models
a <- r.squaredGLMM(L5_MarAprRunoff_MarAprDegree_Ad)
b <- r.squaredGLMM(L5_MarAprRunoff_MarAprDegree_Int)
c <- r.squaredGLMM(L5_MarAprDegree_only)
d <- r.squaredGLMM(L5_MarAprRunoff_only)

a2 <- r.squaredGLMM(L95_MarJunRunoff_MarAprDegree_Ad)
b2 <- r.squaredGLMM(L95_MarJunRunoff_MarAprDegree_Int)
c2 <- r.squaredGLMM(L95_MarJunRunoff_only)
d2 <- r.squaredGLMM(L95_MarAprDegree_only)

a3 <- r.squaredGLMM(L_MarJunRunoff_MarAprDegree_Ad)
b3 <- r.squaredGLMM(L_MarJunRunoff_MarAprDegree_Int)
c3 <- r.squaredGLMM(L_MarJunRunoff_only)
d3 <- r.squaredGLMM(L_MarAprDegree_only)


r2 <- as.data.frame(rbind(a, b, c, d, a2, b2, c2, d2, a3, b2, c3, d3)) %>%
  cbind(data.frame("Model" = c('L5_MarAprRunoff_MarAprDegree_Ad', 'L5_MarAprRunoff_MarAprDegree_Int', 'L5_MarAprDegree_only', 'L5_MarAprRunoff_only', 'L95_MarJunRunoff_MarAprDegree_Ad', 'L95_MarJunRunoff_MarAprDegree_Int', 'L95_MarJunRunoff_only', 'L95_MarAprDegree_only', 'L_MarJunRunoff_MarAprDegree_Ad', 'L_MarJunRunoff_MarAprDegree_Int', 'L_MarJunRunoff_only', 'L_MarAprDegree_only')))

AIC_table_w_R2 <- left_join(AIC_table, r2, by = 'Model')

write.csv(AIC_table_w_R2, '../output/Tables_appendix/Kastl_2022_AIC_table.csv')
```




MY PARTIAL DEPENDENCY PLOTS
```{r}
#library(partial.plot)
library(sjPlot)
library(glmmTMB)

# PLOT MARGINAL EFFECTS
a <- get_model_data(L5_MarAprRunoff_MarAprDegree_Ad, type = "pred", terms = "Degree_Days_sum_scaled_MarApr")
a$Degree_Days_sum_scaled_MarApr <- a$x * 100
a$predicted <- as.Date(a$predicted, origin = "1970-01-01 00:00:00")
a$conf.high <- as.Date(a$conf.high, origin = "1970-01-01 00:00:00")
a$conf.low <- as.Date(a$conf.low, origin = "1970-01-01 00:00:00")

a1 <- ggplot() +
  geom_line(data=a, aes(x=Degree_Days_sum_scaled_MarApr, y=predicted), color='black') + 
  geom_rug(data=Output_StYr_all_lmer_MarApr, aes(x=Degree_Days_sum_scaled_MarApr*100, y=as.Date(Five_pct, origin = "1970-01-01 00:00:00")), alpha = .3, position = "jitter") +
  geom_ribbon(data=a, aes(x=Degree_Days_sum_scaled_MarApr, ymin= a$conf.low,ymax= a$conf.high, stat="identity", alpha=0.05), fill='tomato') +
  labs(x='Degree-days (°C)', y='') +
  theme_bw() +
  theme(legend.position = "none",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        strip.background = element_rect(fill="white"))


b <- get_model_data(L5_MarAprRunoff_MarAprDegree_Ad, type = "pred", terms = "MarApr_Trib_Runoff_mmPerDay")
b$MarApr_Trib_Runoff_mmPerDay <- b$x
b$predicted <- as.Date(b$predicted, origin = "1970-01-01 00:00:00")
b$conf.high <- as.Date(b$conf.high, origin = "1970-01-01 00:00:00")
b$conf.low <- as.Date(b$conf.low, origin = "1970-01-01 00:00:00")

b1 <- ggplot() +
  geom_line(data=b, aes(x=MarApr_Trib_Runoff_mmPerDay, y=predicted), color='black') + 
  geom_rug(data=Output_StYr_all_lmer_MarApr, aes(x=MarApr_Trib_Runoff_mmPerDay, y=as.Date(Five_pct, origin = "1970-01-01 00:00:00")), alpha = .3, position = "jitter") +
  geom_ribbon(data=b, aes(x=MarApr_Trib_Runoff_mmPerDay, ymin= b$conf.low,ymax= b$conf.high, stat="identity", alpha=0.05), fill='steel blue') +
  labs(y='Outmigration start date', x='Daily runoff* (mm)') +
  theme_bw() +
  theme(legend.position = "none",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        strip.background = element_rect(fill="white"))


c <- get_model_data(L95_MarJunRunoff_MarAprDegree_Ad, type = "pred", terms = "Degree_Days_sum_scaled_MarApr")
c$Degree_Days_sum_scaled_MarApr <- c$x * 100
c$predicted <- as.Date(c$predicted, origin = "1970-01-01 00:00:00")
c$conf.high <- as.Date(c$conf.high, origin = "1970-01-01 00:00:00")
c$conf.low <- as.Date(c$conf.low, origin = "1970-01-01 00:00:00")

c1 <- ggplot() +
  geom_line(data=c, aes(x=Degree_Days_sum_scaled_MarApr, y=predicted), color='black') + 
  geom_rug(data=Output_StYr_all_lmer_MarApr, aes(x=Degree_Days_sum_scaled_MarApr*100, y=as.Date(NinetyFive_pct, origin = "1970-01-01 00:00:00")), alpha = .3, position = "jitter") +
  geom_ribbon(data=c, aes(x=Degree_Days_sum_scaled_MarApr, ymin= c$conf.low,ymax= c$conf.high, stat="identity", alpha=0.05), fill='tomato') +
  labs(x='Degree-days (°C)', y='') +
  theme_bw() +
  theme(legend.position = "none",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        strip.background = element_rect(fill="white"))


d <- get_model_data(L95_MarJunRunoff_MarAprDegree_Ad, type = "pred", terms = "MarJun_Trib_Runoff_mmPerDay")
d$MarJun_Trib_Runoff_mmPerDay <- d$x
d$predicted <- as.Date(d$predicted, origin = "1970-01-01 00:00:00")
d$conf.high <- as.Date(d$conf.high, origin = "1970-01-01 00:00:00")
d$conf.low <- as.Date(d$conf.low, origin = "1970-01-01 00:00:00")

d1 <- ggplot() +
  geom_line(data=d, aes(x=MarJun_Trib_Runoff_mmPerDay, y=predicted), color='black') + 
  geom_rug(data=Output_StYr_all_lmer_MarApr, aes(x=MarJun_Trib_Runoff_mmPerDay, y=as.Date(NinetyFive_pct, origin = "1970-01-01 00:00:00")), alpha = .3, position = "jitter") +
  geom_ribbon(data=d, aes(x=MarJun_Trib_Runoff_mmPerDay, ymin= d$conf.low,ymax= d$conf.high, stat="identity", alpha=0.05), fill='steel blue') +
  labs(y='Outmigration end date', x='Daily runoff^ (mm)') +
  theme_bw() +
  theme(legend.position = "none",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        strip.background = element_rect(fill="white"))


e <- get_model_data(L_MarJunRunoff_MarAprDegree_Ad, type = "pred", terms = "MarJun_Trib_Runoff_mmPerDay")
e$MarJun_Trib_Runoff_mmPerDay <- e$x

e1 <- ggplot() +
  geom_line(data=e, aes(x=MarJun_Trib_Runoff_mmPerDay, y=predicted), color='black') + 
  geom_rug(data=Output_StYr_all_lmer_MarApr, aes(x=MarJun_Trib_Runoff_mmPerDay, y=as.Date(Dur_5_95, origin = "1970-01-01 00:00:00")), alpha = .3, position = "jitter") +
  geom_ribbon(data=e, aes(x=MarJun_Trib_Runoff_mmPerDay, ymin= e$conf.low,ymax= e$conf.high, stat="identity", alpha=0.05), fill='steel blue') +
  labs(y='Outmigration window (days)', x='Daily runoff^ (mm)') +
  theme_bw() +
  theme(legend.position = "none",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        strip.background = element_rect(fill="white"))


Compilation_PartialDep <- cowplot::plot_grid(b1, a1, d1, c1, e1, labels = c("(a)","(b)","(c)","(d)","(e)"), ncol = 2)
Compilation_PartialDep
ggsave("../output/Figs_main_paper/Kastl_2022_partial_dependency_plots.jpg",  width=8, height=9, dpi=600)
```





*ADDITIVE MODELS* PLOT CONFIDENCE INTERVAL WHISKER PLOT
```{r}

############ Extract effect sizes and R2 from model ###############

# RUNOFF FIGURE
# fixed effect coefficients
fix_5r <- as.data.frame(fixef(L5_MarAprRunoff_MarAprDegree_Ad))
CI_5r <- as.data.frame(confint(L5_MarAprRunoff_MarAprDegree_Ad, method='Wald')) # get walds confidence interval
#CI2 <- confint(m1,method="boot") # bootstrap, time intensive, similar results to Wald method, see Bolker GLMM worked examples. 
# Effect Sizes from model averaged coefficients
effectsize_5r <- data.frame(variable = rownames(fix_5r),
                         coef = fix_5r[ ,1],
                         ci.low = CI_5r[-c(1:2),1],
                         ci.high = CI_5r[-c(1:2),2]) %>%
  filter(variable == 'MarApr_Trib_Runoff_mmPerDay')

############## Plot effect sizes and relative importance ##################
#effectsize <- effectsize[-1,] # removing intercept effect size
# # # assign new names for variables in order of importance # check order of importance # # #
effectsize_5r$name <- 'Daily Runoff*'
effectsize_5r$signif <- ifelse(effectsize_5r$ci.low<0 & effectsize_5r$ci.high>0, 'no', 'yes')


fix_95r <- as.data.frame(fixef(L95_MarJunRunoff_MarAprDegree_Ad))
CI_95r <- as.data.frame(confint(L95_MarJunRunoff_MarAprDegree_Ad, method='Wald')) # get walds confidence interval
#CI2 <- confint(m1,method="boot") # bootstrap, time intensive, similar results to Wald method, see Bolker GLMM worked examples. 

# Effect Sizes from model averaged coefficients
effectsize_95r <- data.frame(variable = rownames(fix_95r),
                         coef = fix_95r[ ,1],
                         ci.low = CI_95r[-c(1:2),1],
                         ci.high = CI_95r[-c(1:2),2]) %>%
  filter(variable == 'MarJun_Trib_Runoff_mmPerDay')

effectsize_95r$name <- 'Daily Runoff^'
effectsize_95r$signif <- ifelse(effectsize_95r$ci.low<0 & effectsize_95r$ci.high>0, 'no', 'yes')


fix_5_95r <- as.data.frame(fixef(L_MarJunRunoff_MarAprDegree_Ad))
CI_5_95r <- as.data.frame(confint(L_MarJunRunoff_MarAprDegree_Ad, method='Wald')) # get walds confidence interval

# Effect Sizes from model averaged coefficients
effectsize_5_95r <- data.frame(variable = rownames(fix_5_95r),
                         coef = fix_5_95r[ ,1],
                         ci.low = CI_5_95r[-c(1:2),1],
                         ci.high = CI_5_95r[-c(1:2),2]) %>%
  filter(variable == 'MarJun_Trib_Runoff_mmPerDay')

effectsize_5_95r$name <- 'Daily Runoff^'
effectsize_5_95r$signif <- ifelse(effectsize_95r$ci.low<0 & effectsize_5_95r$ci.high>0, 'no', 'yes')



# DEGREE-DAY FIGURE
# fixed effect coefficients
fix_5d <- as.data.frame(fixef(L5_MarAprRunoff_MarAprDegree_Ad))
CI_5d <- as.data.frame(confint(L5_MarAprRunoff_MarAprDegree_Ad, method='Wald')) # get walds confidence interval

# Effect Sizes from model averaged coefficients
effectsize_5d <- data.frame(variable = rownames(fix_5d),
                         coef = fix_5d[ ,1],
                         ci.low = CI_5d[-c(1:2),1],
                         ci.high = CI_5d[-c(1:2),2]) %>%
  filter(variable == 'Degree_Days_sum_scaled_MarApr')

# FOR TABLE (I need actual estimates, unlike the whisker figure that keeps degree-days scaled):
# Divide "Estimate" and "Std Error" by 100 because I scaled Degree_Days_sum_scaled_MarApr, dividing it by 100. 
# Explanation: y = mx + nz + b. m = Estimate. x = degree-days. I initially divided x by 100, so I need to divide m by 100 because  m*(x/100)=x*(m/100) 
effectsize_5d_table <- data.frame(variable = rownames(fix_5d),
                         coef = (fix_5d[ ,1])/100,
                         ci.low = (CI_5d[-c(1:2),1])/100,
                         ci.high = (CI_5d[-c(1:2),2])/100) %>%
  filter(variable == 'Degree_Days_sum_scaled_MarApr')

############## Plot effect sizes and relative importance ##################

# # # assign new names for variables in order of importance # check order of importance # # #
effectsize_5d$name <- 'Degree-days'
effectsize_5d$signif <- ifelse(effectsize_5d$ci.low<0 & effectsize_5d$ci.high>0, 'no', 'yes')

effectsize_5d_table$name <- 'Degree-days'
effectsize_5d_table$signif <- ifelse(effectsize_5d$ci.low<0 & effectsize_5d$ci.high>0, 'no', 'yes')

fix_95d <- as.data.frame(fixef(L95_MarJunRunoff_MarAprDegree_Ad))
CI_95d <- as.data.frame(confint(L95_MarJunRunoff_MarAprDegree_Ad, method='Wald')) # get walds confidence interval

# Effect Sizes from model averaged coefficients
effectsize_95d <- data.frame(variable = rownames(fix_95d),
                         coef = fix_95d[ ,1],
                         ci.low = CI_95d[-c(1:2),1],
                         ci.high = CI_95d[-c(1:2),2]) %>%
  filter(variable == 'Degree_Days_sum_scaled_MarApr')

effectsize_95d$name <- 'Degree-days'
effectsize_95d$signif <- ifelse(effectsize_95d$ci.low<0 & effectsize_95d$ci.high>0, 'no', 'yes')



# FOR TABLE (I need actual estimates, by scaling back, unlike the whisker figure that keeps degree-days scaled):
# Divide "Estimate" and "Std Error" by 100 because I scaled Degree_Days_sum_scaled_MarApr, dividing it by 100. 
# Explanation: y = mx + nz + b. m = Estimate. x = degree-days. I initially divided x by 100, so I need to divide m by 100 because m*(x/100)=x*(m/100)

effectsize_95d_table <- data.frame(variable = rownames(fix_95d),
                         coef = (fix_95d[ ,1])/100,
                         ci.low = (CI_95d[-c(1:2),1])/100,
                         ci.high = (CI_95d[-c(1:2),2])/100) %>%
  filter(variable == 'Degree_Days_sum_scaled_MarApr')

effectsize_95d_table$name <- 'Degree-days'
effectsize_95d_table$signif <- ifelse(effectsize_95d$ci.low<0 & effectsize_95d$ci.high>0, 'no', 'yes')





fix_5_95d <- as.data.frame(fixef(L_MarJunRunoff_MarAprDegree_Ad))
CI_5_95d <- as.data.frame(confint(L_MarJunRunoff_MarAprDegree_Ad, method='Wald')) # get walds confidence interval

# Effect Sizes from model averaged coefficients
effectsize_5_95d <- data.frame(variable = rownames(fix_5_95d),
                         coef = fix_5_95d[ ,1],
                         ci.low = CI_5_95d[-c(1:2),1],
                         ci.high = CI_5_95d[-c(1:2),2]) %>%
  filter(variable == 'Degree_Days_sum_scaled_MarApr')

effectsize_5_95d$name <- 'Degree-days^'
effectsize_5_95d$signif <- ifelse(effectsize_95d$ci.low<0 & effectsize_5_95d$ci.high>0, 'no', 'yes')





effectsize_all <- rbind(effectsize_5d, effectsize_5r, effectsize_95r, effectsize_95d, effectsize_5_95r, effectsize_5_95d)
effectsize_all_table <- rbind(effectsize_5d_table, effectsize_5r, effectsize_95r, effectsize_95d_table, effectsize_5_95r, effectsize_5_95d)


p1 <- ggplot(effectsize_all, aes(coef, name, colour=signif)) 
peffect <- p1 + 
  geom_vline(xintercept = 0,linetype='dashed') +
  geom_point(size =3, stat = "identity") + 
  scale_color_manual(values=c('black', 'dark grey')) +
  geom_errorbarh(aes(xmin=ci.low, xmax=ci.high), height=0.3, size=1) +
  scale_x_continuous(limits=c(-25,25), name='Coefficient (95% C.I.)') +
  theme_bw() +
  theme(axis.title.y=element_blank(),
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        strip.background = element_rect(fill="white"))


# SAME AS ABOVE FIGURE, BUT SEPARATE INTO SECTIONS BY MODEL 
effectsize5dr <- rbind(effectsize_5d, effectsize_5r)
effectsize95dr <- rbind(effectsize_95d, effectsize_95r)
effectsize_5_95dr <- rbind(effectsize_5_95d, effectsize_5_95r)

p2 <- ggplot(effectsize5dr, aes(coef, name, colour=name)) 
peffect2 <- p2 + 
  geom_vline(xintercept = 0,linetype='dashed') +
  geom_point(size =3, stat = "identity") + 
  scale_color_manual(values=c('steel blue', 'tomato')) +
  geom_errorbarh(aes(xmin=ci.low, xmax=ci.high), height=0.2, size=0.5) +
  scale_x_continuous(limits=c(-25,25), name='Coefficient (95% C.I.)') +
  theme_bw() +
  theme(axis.title.y=element_blank(),
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        strip.background = element_rect(fill="white")) +
  labs(subtitle="Outmigration start date") +
  theme(plot.subtitle=element_text(size=10, vjust=-1)) # +
  #annotate(geom='text', label = expression(paste("    Conditional ", R^2 ,": 0.67")), x=-Inf, y=Inf, hjust=0, vjust=1.5, size=2.5)
peffect2

p3 <- ggplot(effectsize95dr, aes(coef, name, colour=name)) 
peffect3 <- p3 + 
  geom_vline(xintercept = 0,linetype='dashed') +
  geom_point(size =3, stat = "identity") + 
  scale_color_manual(values=c('steel blue', 'tomato')) +
  geom_errorbarh(aes(xmin=ci.low, xmax=ci.high), height=0.2, size=0.5) +
  scale_x_continuous(limits=c(-25,25), name='Coefficient (95% C.I.)') +
  theme_bw() +
  theme(axis.title.y=element_blank(),
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        strip.background = element_rect(fill="white")) +
    labs(subtitle="Outmigration end date") +
    theme(plot.subtitle=element_text(size=10, vjust=-1)) # +
  #annotate(geom='text', label = expression(paste("    Conditional ", R^2 ,": 0.67")), x=-Inf, y=Inf, hjust=0, vjust=1.5, size=2.5)
peffect3

p4 <- ggplot(effectsize_5_95dr, aes(coef, name, colour=name)) 
peffect4 <- p4 + 
  geom_vline(xintercept = 0,linetype='dashed') +
  geom_point(size =3, stat = "identity") + 
  scale_color_manual(values=c('steel blue', 'tomato')) +
  geom_errorbarh(aes(xmin=ci.low, xmax=ci.high), height=0.2, size=0.5) +
  scale_x_continuous(limits=c(-25,25), name='Coefficient (95% C.I.)') +
  theme_bw() +
  theme(axis.title.y=element_blank(),
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        strip.background = element_rect(fill="white")) +
    labs(subtitle="Outmigration window") +
    theme(plot.subtitle=element_text(size=10, vjust=-1)) #+
  #annotate(geom='text', label = expression(paste("    Conditional ", R^2 ,": 0.48")), x=-Inf, y=Inf, hjust=0, vjust=1.5, size=2.5)
peffect4


Composite_Coef <- cowplot::plot_grid(peffect2 + theme(axis.text.x = element_blank(),
                                    axis.ticks.x = element_blank(),
                                    axis.title.x = element_blank() ),
                                    peffect3 + 
                                      theme(axis.text.x = element_blank(),
                                    axis.ticks.x = element_blank(),
                                    axis.title.x = element_blank() ),
                                    peffect4 +
                                      theme(axis.title.x = element_text(size=10) ),
                                    nrow=3,
                                    labels = c("(a)","(b)","(c)"))
Composite_Coef
ggsave("../output/Figs_main_paper/Kastl_2022_coefficients_effect4.jpg",  width=5, height=5, dpi=600)
```


Calculate the change in mean outmigration end date, based on the expected water temperature warming by the 2050's (Manuscript Discussion)
```{r}
coef_95_Deg <- effectsize95dr %>% 
  filter(variable == 'Degree_Days_sum_scaled_MarApr') %>%
  select(coef)

# Change in number of days = degrees C expected water temperature increase * 61 days (March-April degree-days) * modeling coefficient * 
new_end <- 1.9 * 61 * (coef_95_Deg/100)

# Add number of seconds of new_end to mean_end
new_end_date <- as.POSIXct(mean_end) + as.numeric(new_end*86400)
```


Calculating the range of predictor variable output over the range of conditions
```{r}
# Duration low-flow prediction
(fix_5_95r['MarJun_Trib_Runoff_mmPerDay',1]*0.1144884) + fix_5_95r['(Intercept)',1]

# Duration low-flow CI range
(CI_5_95r['MarJun_Trib_Runoff_mmPerDay','2.5 %']*0.1144884) + CI_5_95r['(Intercept)','2.5 %']
(CI_5_95r['MarJun_Trib_Runoff_mmPerDay','97.5 %']*0.1144884) + CI_5_95r['(Intercept)','97.5 %']

#Duration high-flow prediction
(fix_5_95r['MarJun_Trib_Runoff_mmPerDay',1]*1.7602132) + fix_5_95r['(Intercept)',1]

# Duration high-flow CI range
(CI_5_95r['MarJun_Trib_Runoff_mmPerDay','2.5 %']*1.7602132) + CI_5_95r['(Intercept)','2.5 %']
(CI_5_95r['MarJun_Trib_Runoff_mmPerDay','97.5 %']*1.7602132) + CI_5_95r['(Intercept)','97.5 %']  
```


Show p-values for top models, within ~2.0 of one another (check the lowest AIC values from above for top models)
Lower (in this case more negative) AIC values indicates superior model fit.
```{r}
library(MuMIn)
library(gdata)

# FOR THE INTERACTION TERM MODELS
summary(L5_MarAprRunoff_MarAprDegree_Int)$coefficients

summary(L95_MarJunRunoff_MarAprDegree_Int)$coefficients






# FOR ORIGINAL ADDITIVE MODELS
summary(L5_MarAprRunoff_MarAprDegree_Ad)$coefficients

summary(L95_MarJunRunoff_MarAprDegree_Ad)$coefficients

summary(L_MarJunRunoff_only)$coefficients
summary(L_MarJunRunoff_MarAprDegree_Ad)$coefficients


# Divide "Estimate" and "Std Error" by 100 because I scaled Degree_Days_sum_scaled_MarApr, dividing it by 100. 
# Explanation: y = mx + nz + b. m = Estimate. x = degree-days. I initially divided x by 100, so I need to divide m by 100 because m*(x/100)=x*(m/100) 
Estimate_L5_MarAprRunoff_MarAprDegree_Ad <- (summary(L5_MarAprRunoff_MarAprDegree_Ad)$coefficients['Degree_Days_sum_scaled_MarApr', 'Estimate'])/100
Estimate_L5_MarAprRunoff_MarAprDegree_Ad

StdError_L5_MarAprRunoff_MarAprDegree_Ad <- (summary(L5_MarAprRunoff_MarAprDegree_Ad)$coefficients['Degree_Days_sum_scaled_MarApr', 'Std. Error'])/100
StdError_L5_MarAprRunoff_MarAprDegree_Ad


Estimate_L95_MarAprDegree_only <- (summary(L95_MarAprDegree_only)$coefficients['Degree_Days_sum_scaled_MarApr', 'Estimate'])/100
Estimate_L95_MarAprDegree_only

StdError_L95_MarAprDegree_only <- (summary(L95_MarAprDegree_only)$coefficients['Degree_Days_sum_scaled_MarApr', 'Std. Error'])/100
StdError_L95_MarAprDegree_only


Estimate_L_MarJunRunoff_MarAprDegree_Ad <- (summary(L_MarJunRunoff_MarAprDegree_Ad)$coefficients['Degree_Days_sum_scaled_MarApr', 'Estimate'])/100
Estimate_L_MarJunRunoff_MarAprDegree_Ad

StdError_L_MarJunRunoff_MarAprDegree_Ad <- (summary(L_MarJunRunoff_MarAprDegree_Ad)$coefficients['Degree_Days_sum_scaled_MarApr', 'Std. Error'])/100
StdError_L_MarJunRunoff_MarAprDegree_Ad
```

