w_steelhead_July_ISU.Rmd

---
title: Columbia Wild Summer Steelhead Inseason Forecast for Counts 7/1-7/31 @ Bonneville Dam
author: Thomas Buehrens (tbuehrens@dfw.wa.gov) 
output:
  html_document:
    code_folding: hide
    fig_caption: yes
    theme: cerulean
    toc: yes
    toc_depth: 3
    toc_float: yes
  pdf_document:
    toc: yes
    toc_depth: '3'
---
<script>
   $(document).ready(function() {
     $head = $('#header');
     $head.prepend('<img src=\"https://privatelands.wdfw.wa.gov/wdfwlogo_clrnotxt.png"\" style=\"float: right;width: 150px;\"/>')
   });
</script>

***

Last Updated `r format(Sys.time(), '%m/%d/%Y')`.

***

# Overview
This script fits Seasonal Auto-regressive Integrated Moving Average models to in-season salmon return data to predict total return. It does this by aggregating the run-size into two periods: counts prior to a particular date (often the date of the most recent observation), and counts after that date (i.e., the remainder of the run). It then predicts the remainder of the run using the pattern of abundance before and after that date in previous years, in conjunction with covariates. This software evaluates the model's skill at predicting the remainder of the runsize based on data collected up to a particular date. A forecast for the final year, including confidence intervals, is produced.

<!-- ## Setup -->
<!-- All analyses require R software [**(link)**](https://cran.r-project.org/) (v3.4.3) for data retrieval, data processing, and summarizing model results. Here we configure R to perform our analysis and generate our outputs -->
```{r set_options, echo = TRUE, message = FALSE}
options(width = 100)
knitr::opts_chunk$set(message = FALSE)
set.seed(123)
```

<!-- We also need a couple of helper functions which we will define -->
```{r load_funcs, message = FALSE, warning = FALSE,results = "hide"}
wd_functions<-"functions"
sapply(FUN = source, paste(wd_functions, list.files(wd_functions), sep="/"))
```

<!-- Here we will load & install packages we need to use (needs internet connection if packages not already installed) -->
```{r load_packages, message = FALSE, warning = FALSE,results = "hide"}
packages_list<-c("tidyverse"
                 ,"forecast"
                 ,"mgcv"
                 ,"ggplot2"
                 ,"MASS"
                 ,"RColorBrewer"
                 ,"kableExtra"
                 ,"lubridate"
                 ,"modelr"
                 ,"kableExtra"
                 ,"reshape2"
                 ,"ggfortify"
                 ,"clock"
                 ,"smooth"
                 ,"scales"
                 ,"gtools"
                 ,"CombMSC"
                 ,"here"
                 ,"MuMIn"
                 )
install_or_load_pack(pack = packages_list)
```

## User Inputs 
Look at this section to see user inputs regarding the years and dates of data analyzed and covariates used.
```{r user_inputs_data, message = FALSE, warning = FALSE,results = "show"}
#=========
# Raw Data
#=========
yr_start<-2003
yr_end<-year(Sys.Date())
dam="BON"
species="WStlhd"
#===========================
# Summarization for analysis
#===========================
date_start_analysis<-ymd("2003/6/01") #this will be the first date of data used in the analysis (and this month/day first in years thereafter)
date_end_analysis<-ymd("2022/7/31") #this will be the last date of data used in the analysis (and this month/day last in years preceding)
forecast_period_start_m<-7 #this will be the month associated with the first month/day of the seasonal estimate period each year...can be after first data used to estimate it or the same
forecast_period_start_d<-1 #this will be the month day associated with the first month/day of the seasonal estimate period each year...can be after first data used to estimate it or the same
use_freshest_data = T #use all data up to "today" or only up to beginning of forecast period
covariates<-c("lag1_log_SAR1","lag2_NPGO", "lag1_log_SAR2","lag1_NPGO")
#==================
#forecasting params
#==================
leave_yrs<- 15
covariates<-c("lag1_log_SAR1","lag1_log_SAR2","lag1_NPGO","lag2_NPGO")
p1_covariates_only =c("zl_flow")
plot_results = F
first_forecast_period = 2
write_model_summaries = TRUE
find_best=T
stack_metric = "MAPE"
```

## Get Raw Data
Data is loaded here a a snapshot of the data used in the analysis (before aggregation into the two periods described) is provided.
```{r get_data, message=FALSE, warning=FALSE, results="show"}

dat<-get_dam_data(
  yr_start = yr_start,
  yr_end = yr_end,
  dam = dam,
  species = species
)

#=========================================================
#get NPGO data
#=========================================================
NPGO<-read_table("http://www.o3d.org/npgo/npgo.php",skip=29,col_names=F,comment="#")%>%
  filter(!is.na(X2))%>%
  dplyr::rename(Year=X1,Month=X2,NPGO=X3)%>%
  mutate(Year=as.numeric(Year))%>%
  group_by(Year)%>%
  add_tally()%>%
  #filter(!Month>6)%>% #use only spring (Jan-June) NPGO
  filter(!n < 12)%>% #use only complete years
  group_by(Year)%>%
  dplyr::summarise(NPGO=mean(NPGO))%>%
  mutate(Year=Year+1, lag1_NPGO = NPGO,lag2_NPGO = lag(NPGO))%>%
  dplyr::select(year=Year,lag1_NPGO,lag2_NPGO)

# #=========================================================
# #get PIT tag survival/age data
# #=========================================================
PIT<-read_csv("https://raw.githubusercontent.com/tbuehrens/Salmon_Forecast_Example/main/data/Columbia_Steelhead_SAR_DART_11_1_2021.csv")%>%
  dplyr::rename(OutmigrationYear=year)%>%
  mutate(Year=OutmigrationYear+2)%>%
  filter(Pop=="SNA")

PIT<-PIT%>%bind_cols(data.frame(SAR1=gam(cbind(ocean1Count,juvCount-ocean1Count)~s(OutmigrationYear,k=(dim(PIT)[1])),family=binomial,data=PIT)$fitted))%>%
  bind_cols(data.frame(SAR2=c(gam(cbind(ocean2Count,juvCount-ocean2Count)~s(OutmigrationYear,k=(dim(PIT)[1]-1)),family=binomial,data=PIT)$fitted,NA)))%>%
  mutate(lag1_log_SAR1 = log(SAR1),lag1_log_SAR2=lag(log(SAR2),1))%>%
  dplyr::select(year=Year,lag1_log_SAR1,lag1_log_SAR2)

dat<-dat%>%
  left_join(NPGO)%>%
  left_join(PIT)

print(head(dat))
```

<!-- ## Summarize Data for Analysis -->
```{r summarize_data, message=FALSE, warning=FALSE, results="show"}
summarized_data<-prepare_data(series = dat,
                  date_start_analysis = date_start_analysis,
                  date_end_analysis = date_end_analysis,
                  forecast_period_start_m = forecast_period_start_m, #inclusive 
                  forecast_period_start_d = forecast_period_start_d, #inclusive
                  use_freshest_data = use_freshest_data,
                  covariates = covariates,
                  p1_covariates_only=p1_covariates_only
                  )

```

## Results
This section present a a results table showing the performance of the prediction model in previous years as well as a forecast in the current year. The results in the table are graphed below.
```{r Analysis_v2, message=FALSE, warning=FALSE, results="show"}
if(find_best ==T){
  best_covariates<-all_subsets(series=summarized_data$series,covariates=covariates,min=0,max=4)
  saveRDS(best_covariates,"best_covariates_steelhead_july.rds")
}


# results<-inseason_forecast(series=summarized_data$series,
#                   leave_yrs=leave_yrs,
#                   covariates= c("lag2_NPGO","lag1_log_SAR2","lag2_log_smolts","var_flow","zl_flow"), #use to automate variable selection: best_covariates[[1]][[best_covariates[[2]]$model_num[1]]],
#                   first_forecast_period = first_forecast_period,
#                   plot_results = plot_results,
#                   write_model_summaries = write_model_summaries,
#                   forecast_period_start_m =  forecast_period_start_m, #inclusive 
#                   forecast_period_start_d =  forecast_period_start_d, #inclusive
#                   obs_period_2 = summarized_data$obs_period_2,
#                   p1_covariates_only = p1_covariates_only
#                   )

best_covariates<-readRDS("best_covariates_steelhead_july.rds")

best_covariates[[2]]%>%
  slice(1:10)%>%
  kbl(caption = "Table 1. Covariate Model Selection Results.",digits =3)%>%
  kable_classic(full_width = F, html_font = "Cambria")

results<-inseason_forecast_v2(series=summarized_data$series,
                  leave_yrs=leave_yrs,
                  covariates= best_covariates[[1]][best_covariates[[2]]$model_num[1:10]],
                  first_forecast_period = first_forecast_period,
                  plot_results = plot_results,
                  write_model_summaries = write_model_summaries,
                  forecast_period_start_m =  forecast_period_start_m, #inclusive 
                  forecast_period_start_d =  forecast_period_start_d, #inclusive
                  obs_period_2 = summarized_data$obs_period_2,
                  p1_covariates_only = p1_covariates_only,
                  stack_metric = stack_metric
                  )

model_list<-lapply(best_covariates[[1]][best_covariates[[2]]$model_num[1:10]],function(x) paste(x,collapse = " + "))%>%
  unlist()%>%
  as_tibble()%>%
  add_rownames()%>%
  dplyr::rename(model=rowname,model_name=value)

results%>%
  filter(year==max(year))%>%
  left_join(model_list)%>%
  mutate(model_name = ifelse(is.na(Stacking_weight),"ensemble",model_name))%>%
  dplyr::select(year,model_name,predicted_abundance,`Lo 50`,`Hi 50`,`Lo 95`,`Hi 95`,Stacking_weight)%>%
  kbl(caption = "Table: Individual Ensemble Member Predictions and Weights",digits =2)%>%
  kable_classic(full_width = F, html_font = "Cambria")


if(length(unique(results$model))>1){

  results_best<-results%>%
  ungroup()%>%
  filter(model =="ensemble")
}


p_2_start_m<-summarized_data$p_2_start_m
p_2_start_d<-summarized_data$p_2_start_d


results_best%>%
  dplyr::select(!c("period","train_test"))%>%
  kbl(caption = paste0("Table 2. Forecasts of ",dam," dam counts of ",species, " from ",forecast_period_start_m,"/",forecast_period_start_d,"-",month(date_end_analysis),"/",mday(date_end_analysis)," using only past years' data for the forecast period and in-season counts from ", month(date_start_analysis),"/",mday(date_start_analysis)," through ", month(max(dat$date)),"/",mday(max(dat$date)), " in each year. Additional covariates are included in the table below."),digits =2)%>%
  kable_classic(full_width = F, html_font = "Cambria")

print(paste0("Mean Absolute Percent Error (MAPE) = ",mean(abs((results_best$error/results$abundance)*100),na.rm = T)))
print(paste0("Median Symmetric Accuracy; MSA) = ",100*(exp(median(abs(log(results_best$predicted_abundance/results$abundance)),na.rm = T))-1)))

p<-ggplot(results_best,aes(x=year,y=predicted_abundance))+
  geom_ribbon(aes(ymin=`Lo 95`,ymax=`Hi 95`),color=NA,alpha=0.5,fill = "cadetblue")+
  geom_ribbon(aes(ymin=`Lo 50`,ymax=`Hi 50`),color=NA,alpha=0.5,fill = "cadetblue")+
  geom_line()+
  geom_point(aes(x=year,y=abundance))+
  ylim(0,NA)+
  scale_x_continuous(breaks=unique(results$year))

print(p)
```