This report is automatically generated with the R package knitr (version 1.40) .

---
title: |
       | TLG Eval. (1998-2019 data) & Black Bass Hg Imp. Goal (2000-2019 data) 
       | Test Grouping Data by Subarea vs. Subarea & Year
author: "Mercury Program and Basin Planning Unit"
date: "6/30/2021"
output:
  html_document:
    code_folding: show
    toc: TRUE
    toc_float: TRUE
    toc_depth: 3
runtime: shiny
assets:
  css:
    - "http://fonts.googleapis.com/css?family=Raleway:300"
    - "http://fonts.googleapis.com/css?family=Oxygen"
---

<style>
body{
  font-family: 'Oxygen', sans-serif;
  font-size: 16px;
  line-height: 24px;
}

h1,h2,h3,h4 {
  font-family: 'Raleway', sans-serif;
}

.container { width: 1250px; }
h3 {
  background-color: #D4DAEC;
  text-indent: 50px; 
}
h4 {
  text-indent: 75px;
  margin-top: 35px;
  margin-bottom: 5px;
}
</style>

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo=TRUE, warning=FALSE, message=FALSE, fig.width=9.5)
```

This markdown looks at all fish species data from 1998 through 2019 & Black Bass data from 2000 through 2019. It also tests whether pooling data by subarea or by year and subarea provides a better model. Pooling data by subarea and year provides the lowest avg SER of predicted BB safe concentration models (see Regression 2, Alt 1 & Alt 2 results).   


```{r Libraries, include=FALSE, message=FALSE}
library(knitr)      # generate html tables for this report
library(kableExtra) # better formatting of tables
library(DT)         # creates sortable tables
library(shiny)

# Having issue with runtime:shiny resetting the WD of rproj so need to use setwd()
wd <- rstudioapi::getActiveProject()
setwd(wd)

source("R Functions/functions_estimate NDDNQ values.R") # Load first so select() isn't masked from dplyr
source("R Functions/functions_QA data.R")
source("R Functions/function_predictionModels.R")
source("R Functions/function_predictionModelPlot.R")
source("R Functions/function_TLG_CompositeWeighted_AvgHg.R")
source("R Functions/function_TLG_CompositeWeighted_AvgHg_TissueOptns.R")

# Function to round values in a dataframe
signif_df <- function(x, digits) {
  # round numeric columns
  # x = data frame 
  # digits = number of digits to round
  numeric_columns <- sapply(x, mode) == 'numeric'
  x[numeric_columns] <-  signif(x[numeric_columns], digits)
  x
}

round_df <- function(x, digits) {
  # round numeric columns
  # x = data frame 
  # digits = number of digits to round
  numeric_columns <- sapply(x, mode) == 'numeric'
  x[numeric_columns] <-  round(x[numeric_columns], digits)
  x
}
```

## Load Data
```{r}
setwd(wd)

# Load Fish Data & Supporting Files
fishdf <- loadData("Reeval_Impl_Goals_Linkage_Analysis/Data/Fish/6a_FishMaster_noRepeats.xlsx") 

# Convert columns that are character but should be numeric or date
fishdf <- chara_to_NumDate(fishdf) %>% 
  mutate(Year = year(SampleDate))

# Supproting File - TL Group reference
TLgroups <- readxl::read_xlsx("Reeval_Impl_Goals_Linkage_Analysis/Data/Fish/0_SB FishSpecies-TrophicLevel_Final.xlsx", sheet='Species Trophic Level for R', guess_max = 30000) %>%
  mutate(commonname = tolower(CommonName)) %>%  # create column of TL in lowercase so case isn't an issue for matching names
  mutate(TrophicLevel = ifelse(grepl('na', TrophicLevel), NA, TrophicLevel)) %>% 
  select(commonname, TrophicLevel, TL_minSize) 

TLgroups <- chara_to_NumDate(TLgroups)
```



## Fish Data QC
```{r}
fishdf %>%
  unique_factors(SourceID, Subarea, Analyte, Unit, WBT, CommonName, TissueName)
```

### Change X2 Subarea to West Delta
```{r}
fishData1 <- fishdf %>% 
  mutate(Subarea = recode(Subarea, "X2" = "West Delta")) # Change X2 to West Delta because that's what was done in the 2010 TMDL
```

### Exclude 'mg/kg dw' Units
```{r}
# Select mg/Kg ww (wet weight) because it is the primary unit of measurement for Hg in fish. This is also consistent with the 2010 Linkage. “dw” (dry weight) is provided as a secondary unit of measurement and only included for some samples.
fishData2 <- fishData1 %>% 
  filter(Unit == 'mg/Kg ww')
```

### Add TL Group
```{r}
fishData3 <- fishData2 %>% 
  mutate(commonname = tolower(CommonName)) %>% # create column of TL in lowercase so case isn't an issue for matching names
  left_join(., TLgroups, by='commonname') %>%  # adds TrophicLevel & TL_minSize columns
  select(-commonname) %>% 
  mutate(TrophicLevel = ifelse(!is.na(TL_minSize) & TrophicLevel == 4 & TLAvgLength < TL_minSize, 3, TrophicLevel), # reduce TL4 to TL3 if fish length is smaller than listed min length, based on definitions of TL3 and TL4 fish from "0_SB FishSpecies-TrophicLevel_Final.xlsx”
         TrophicLevel = as.character(TrophicLevel))   # convert to character for graphing

nrow(fishData3) == nrow(fishData2) # this needs to be True

fishData3 %>% 
  select(Subarea, SourceID, SourceRow, StationName, SampleDate, CommonName, TrophicLevel, Result, TLAvgLength, WeightAvg) %>% 
  datatable()
```


## Reevalute 2010 TMDL Staff Report Data Selection Rationale for TLG Evaluation

Text quoted from 2010 TMDL Staff Report Section 4.7.1. Justifications for decisions made by Board staff are in the 2021 TMDL Staff Report.

### Year Range Selection
<em>"First, the data were restricted to samples collected between 1998 and 2001, the period with the most comprehensive sampling across the Delta."</em>
```{r fig.height=5}
fishData3 %>% 
  ggplot(aes(y=Year)) +
  geom_bar(aes(fill=TrophicLevel), width=.8) +
  facet_wrap( ~ Subarea, nrow=1, scales="free_x") +
  scale_y_continuous(minor_breaks=seq(min(fishData2$Year), max(fishData2$Year), 1)) +
  theme_light() 

# Year Selection
fishData4 <- fishData3 %>% 
  filter(Year >= 1998 & year(SampleDate) <= 2019)
```

### Remove/Include Migratory Species
<em>"Second, migratory species (salmon, American shad, steelhead, sturgeon, and striped bass) were excluded. These species likely do not reside year round at the locations in the Delta where they were caught and their tissue mercury levels may not show a positive relationship with the mercury levels in resident animals."</em>
```{r fig.height=5}
migratorySpecies <- c('salmon', 'American shad', 'steelhead', 'sturgeon', 'striped bass')

fishData4 %>% 
  mutate(TrophicLevel= ifelse(grepl(paste(migratorySpecies, collapse="|"), CommonName, ignore.case=T), 'Migratory', TrophicLevel)) %>%
  ggplot(aes(y=Year)) +
  labs(title = 'Showing Migratory Species') +
  geom_bar(aes(fill=TrophicLevel), width=.8) +
  facet_wrap( ~ Subarea, nrow=1, scales="free_x") +
  scale_y_continuous(minor_breaks=seq(min(fishData2$Year), max(fishData2$Year), 1)) +
  theme_light() 

# The same migratory species were removed from the updated reevaluation for the same reasons given in the 2010 Staff Report
fishData4 <- fishData4 %>%
  filter(!grepl(paste(migratorySpecies, collapse="|"), CommonName, ignore.case=T))
```

### Size Range Selection
<em>"Third, fish samples with lengths greater than 500 mm were not included</em> [in the TLG evaluation]<em>. Data for fish larger than 500 mm are available for only some subareas. Capping the size at 500 mm allows comparable data for all Delta subareas."</em>
```{r fig.height=5}
# Original Linkage Data
fishData4 %>% 
  filter(Year >= 1998 & Year <= 2001) %>% 
  ggplot(aes(y=TLAvgLength, x=Subarea)) +
  labs(title = '1998 - 2001 Data') +
  geom_jitter(aes(colour=TrophicLevel), height=0, size=3) +
  geom_hline(yintercept=500, color='red') +
  scale_y_continuous(breaks=seq(0, max(fishData2$TLAvgLength), 200)) +
  theme_light() +
  theme(panel.grid.minor.y = element_blank())
  
# New Data
fishData4 %>% 
  filter(Year > 2001) %>% 
  ggplot(aes(y=TLAvgLength, x=Subarea)) +
  labs(title = '2002 - 2019 Data') +
  geom_jitter(aes(colour=TrophicLevel), height=0, size=3) +
  geom_hline(yintercept=500, color='red') + 
  scale_y_continuous(breaks=seq(0, max(fishData2$TLAvgLength), 200)) +
  theme_light() +
  theme(panel.grid.minor.y = element_blank())

# Size range for TLG evaluation is selected in TLG_CompositeWeighted_AvgHg(). Board staff decided to exclude data above 500mm to maintain consistency with provisions and 2010 TMDL Staff Report.
```

### Tissue Type Selection
<em>"Finally, only fish fillet data were used in the human and eagle trophic level food group analysis. Humans typically consume fish fillets, while wildlife species, including eagles, eat whole fish. However, all the data for large fish typically consumed by eagles and other large wildlife species are from fillet samples, making it necessary to use fillet information for these species. Whole fish data were used for the smaller wildlife species receptors."</em>
```{r fig.height=5}
prep1 <- robinNDDNQ(fishData4, ResultQualCode) # need to est. ND/DNQs so weighted.mean() works correctly
prep2 <-  prep1[[1]]$fitted # select first returned distribution fit 

TLG_AvgHg_tissue <- prep2 %>%
  mutate(TrophicLevel = paste0("TL", TrophicLevel)) %>% 
  TLG_CompositeWeighted_AvgHg_TissueOptns(dataframe=., Subarea) %>% 
  arrange(desc(TLG))

TLG_AvgHg_tissue %>% 
  select(TLG, Tissue, Subarea, CompositeWeightedAvg) %>%
  # distinct(TLG, Subarea, .keep_all=T) %>% 
  pivot_wider(names_from=Subarea, values_from=c(CompositeWeightedAvg)) %>% 
  arrange(desc(TLG), Tissue) %>%
  signif_df(., 3) %>% 
  kbl(caption='Reproducing TMDL Staff Report Table 4.7') %>% 
  kable_paper('hover', full_width=F) %>% 
  kable_styling(fixed_thead = T)

# The table above compares using different tissue selection options based on TLG. To be consistent with the 2010 TMDL Staff Report, Board staff only used fillet data for the human and bald eagle trophic level group analyses. The 2010 TMDL Staff Report only included whole body data for the remaining lower trophic level groups. However, Board staff used fillet and whole body data to include as many lower trophic level groups as possible in the analysis.
```

### Subarea Selection
<em>"Of the eight Delta subareas identified in Section 2.2.2 and Figure 2.2, three of the subareas were not included in the trophic level food group evaluation due to inadequate information. No fish were sampled from the Marsh Creek subarea between 1998 and 2001. In addition, small fish were sampled throughout the Yolo Bypass-South subarea between 1998 and 2001, but large fish were sampled only in the southernmost area; hence, the mercury levels in the trophic level food groups are not geospatially comparable. The only fish sampling conducted in the Yolo Bypass-North subarea took place in Greens Lake, which is not considered representative of the entire subarea. In addition, only large TL4 fish were sampled; no small fish were sampled."</em>
```{r}
fishSampLoc <- fishData4 %>%
  mutate(Species = 'All', 
         uniqName = paste(StationName, Year, round(TargetLatitude,4), round(TargetLongitude,4))) %>%  # needs to be the same code used in 3rd Script when creating variable "AqMaster" in mutate(uniqName = ...)
  distinct(uniqName, .keep_all=T) %>%
  select(Subarea, StationName, Species, Year, TargetLatitude, TargetLongitude, CoordSystem, uniqName) %>% 
  arrange(Subarea, StationName)
```


## Finalize Fish Data Prep
### Estimate ND & DNQ Values
```{r}
fishNDDNQs <- robinNDDNQ(fishData4, ResultQualCode)
```

#### Examine Fitted Distribution Graphs
```{r, echo=FALSE}
ui <- fluidPage(
  selectInput('distName', label='Select Distribution', choices=names(fishNDDNQs), selected=1),
  # textOutput('distText'),
  plotOutput('distPlot')
)

server <- function(input, output, session) {
  # output$distText <- renderText({
  #   paste0('fishNDDNQs[[',input$distName,']]$fitmodel')
  # })
  
  output$distPlot <- renderPlot({
  plot(fishNDDNQs[[input$distName]]$fitmodel) # view plot fits
    title(input$distName, outer=T)
})
}

shinyApp(ui, server, options=list(height=500))
```

#### Select Best Distribution Model
```{r}
fishData5 <- fishNDDNQs[[1]]$fitted  #picked distribution based on plot showing best fit to data

# fishData5 <- randomNDDNQ(fishData4, ResultQualCode)
```

### Check for outliers
```{r, fig.height=5}
# (e.g., orders of magnitude maybe issues with unit conversion)
result_time_Plotly(fishData5, groupByCol=TMDL, showMean=F, interactive=T, logscale=T, showLegend=F)
# Data looks ok, no anomalies

# Look at NDs & DNQs
result_time_Plotly(fishData5 %>% filter(Censored == T),
                   groupByCol=TMDL, showMean=F, interactive=T, logscale=T, showLegend=T)
```

### Save QC'd Data
```{r}
fishData <- fishData5

setwd(wd)
save(fishData, file="Reeval_Impl_Goals_Linkage_Analysis/eval1_TLG Eval & BB Impl Goal/Output/fishData.RData")
```


## TLG Evaluation (Updates 2010 TMDL Staff Report Section 4.7)

### Calcualte TLG Composite Weighted Average for each Subarea
```{r}
# Add "TL" to numberic TrophicLevel - needed for TLG_CompositeWeighted_AvgHg()
TLG <- fishData %>% 
  mutate(TrophicLevel = paste0("TL", TrophicLevel))

# Determine TLG by Subarea
 # Gets the same values as if we grouped by Subarea and then calculated Wt. Avg. by Year
 # TLG_AvgHg_subarea data frame has more info about data used to calculate Composite Wt. Avg. than summary table below
TLG_AvgHg_subarea <- TLG_CompositeWeighted_AvgHg(dataframe=TLG, Subarea) %>% 
  arrange(desc(TLG), Predator)
```

#### View Data Summary
```{r, echo=FALSE}
TLG_AvgHg_shiny <- TLG_AvgHg_subarea %>%
      select(-Predator, -TLGtarget) %>% 
      distinct(TLG, Subarea, .keep_all=T) %>% 
      rename_at(.vars = vars(starts_with("Composite")), # Remove 'Composite' prefix
                .funs = funs(sub("Composite", "", .))) %>% 
      signif_df(., 3)

ui <- fluidPage(
  fluidRow(
    column(width=3, 
           selectInput('subareaName', label='Select Subarea', choices=sort(unique(TLG_AvgHg_shiny$Subarea), decreasing=F), selected=1)
           )
  ),
  tableOutput('subareaTable')
)

server <- function(input, output, session) {

  output$subareaTable <- function() {
    req(input$subareaName)
    
    TLG_AvgHg_shiny %>%
      filter(Subarea == input$subareaName) %>% 
      select(-Subarea) %>% 
      kbl(align='lcccl') %>% 
      kable_paper('hover', full_width=T)
    }
}

shinyApp(ui, server, options=list(height=250))
```

#### Reproduce Table 4.7
Reproduce 2010 TMDL Staff Report Table 4.7 (Mercury Concentrations in Trophic Level Food Groups Sampled in the Delta)
```{r}
Table_4.7 <- TLG_AvgHg_subarea %>% 
  select(TLG, Subarea, CompositeWeightedAvg) %>%
  distinct(TLG, Subarea, .keep_all=T) %>% 
  pivot_wider(names_from=Subarea, values_from=c(CompositeWeightedAvg)) %>% 
  arrange(desc(TLG)) %>% 
  rename(`Trophic Level Group` = TLG)

writexl::write_xlsx(Table_4.7,
                    paste0(wd, '/Reeval_Impl_Goals_Linkage_Analysis/eval1_TLG Eval & BB Impl Goal/Output/TLG Appdx Table X.X 1998-2019 TLG WtAvgs.xlsx'))
  
Table_4.7 %>% 
  signif_df(., 3) %>% 
  kbl(caption='Reproducing TMDL Staff Report Table 4.7') %>% 
  kable_paper('hover', full_width=F) %>% 
  kable_styling(fixed_thead = T)
```


### Join TL4 150-500mm Results to Other TLG Results 
```{r}
# Extract TL4 150-500mm Composite Weighted Avgs for each Subarea
dfTL4_150_500 <- TLG_AvgHg_subarea %>%
  filter(TLG == 'TL4 Fish (150-500 mm)') %>% 
  select(Subarea, CompositeWeightedAvg) %>% 
  distinct(Subarea, .keep_all=T) %>% 
  rename(TL4_150_500 = CompositeWeightedAvg)

# Pair TL4 150-500mm Composite Weighted Avgs with other TLG results
TLG_TL4_150_500 <- TLG_AvgHg_subarea %>% 
  filter(TLG != 'TL4 Fish (150-500 mm)') %>% 
  left_join(.,  dfTL4_150_500, by='Subarea')  # adds TL4_150_500 column of Composite Weighted Avg for each Subarea
```


### Regression: TL4 150-500mm Predicted Concentrations
```{r}
PredictTL4_150_500_models <- TLG_TL4_150_500 %>%  
  group_by(TLG, Predator) %>%
  group_modify(~ predictionModels(.x, .knownYcol=TL4_150_500, .knownXcol=CompositeWeightedAvg, knownX=.x$TLGtarget[1]), preferExtrapModel=NULL) %>%
  rename(PredictTL4_150_500 = predictedY) %>% 
  ungroup
```

#### Regressions: TLG Composite Wt. Avg.
<b>`r paste(nrow(PredictTL4_150_500_models), 'Total Regressions Generated ~ distinct(TLG) would select', nrow(PredictTL4_150_500_models %>% distinct(TLG)), 'regressions (1 for each TLG)' )`</b>
```{r}
# Manually Check Potentially Suspect Regressions (identified by filter() conditions)
PredictTL4_150_500_models%>%
  select(TLG, Predator, ModelName, PredictTL4_150_500, PredictionType, n) %>%
  distinct(TLG, PredictionType, .keep_all=T) %>% 
  filter(grepl("gam", ModelName, ignore.case=T) |
         PredictionType == "Extrapolated") %>% 
  arrange(desc(TLG)) %>% 
  kbl(align='llcc', caption=paste('Check These', nrow(.), 'Suspect Regressions (identified by filter() conditions)')) %>% 
  kable_paper('hover', full_width=T) %>% 
  kable_styling(fixed_thead = T)
```
<em>Model with lowest SER for each TLG & Pred Group is displayed first. Board staff manually checked that the selected models are appropriate. Board staff also checked if keeping extrapolated value is appropriate.</em><br><br>
```{r, echo=FALSE}
ui <- fluidPage(
  fluidRow(
    column(width=3, 
           selectInput('tlgName', label='Select TLG on X-axis', choices=sort(unique(PredictTL4_150_500_models$TLG), decreasing=T), selected=1)
           ),
  column(width=3,
         selectInput('predName', label='Select Pred Group', choices=NULL)
         ),
  column(width=3,
         selectInput('model', label='Select Model', choices=NULL)
         )
  ),
  plotlyOutput('tlgPlot')
)

server <- function(input, output, session) {
  data1 <- reactive({
    PredictTL4_150_500_models %>%
      filter(TLG == input$tlgName)
    })
  observeEvent(data1(), {
    updateSelectInput(session, 'predName',
                      choices=data1()$Predator, selected=data1()$Predator[1])
    })
  
  data2 <- reactive({
    req(input$predName)
    data1() %>%
      filter(Predator == input$predName)
    })
   observeEvent(data2(), {
    updateSelectInput(session, 'model',
                      choices=data2()$ModelName, selected=data2()$ModelName[1])
    })
   
   data3 <- reactive({
    req(input$model)
    data2() %>%
      filter(ModelName == input$model)
   })
    

  output$tlgPlot <- renderPlotly({
    req(input$model)
    
    title  <- paste0('TL4 Fish (150-500 mm) vs. ', input$tlgName)
    xTitle <- input$tlgName
    yTitle <- 'TL4 Fish (150-500mm)'
    knownXlabel <- paste0("TLG Target<sub>",data2()$Predator,"</sub>")
    
    predictionModelPlot(data3(), knownXlabel=knownXlabel,
                        .predictedYcol=PredictTL4_150_500, .groupByCol=Subarea, allGroups=sort(unique(fishData$Subarea)),
                        title=title, xTitle=xTitle, yTitle=yTitle, showLegend=T)
    })
}

shinyApp(ui, server, options=list(height=500))
```

#### Select Models for Predicted TL4 150-500mm
```{r}
handSelectedModels <- PredictTL4_150_500_models %>%
  filter(TLG == 'TL3 Fish (50 - <150 mm)' & ModelName == 'log')

# Replace select lowest SER models with hand selected models
PredictTL4_150_500 <- PredictTL4_150_500_models %>%
  distinct(TLG, Predator, .keep_all=T) %>% 
  anti_join(., handSelectedModels, by=c('TLG')) %>% # remove models that will be replaced by hand selection
  bind_rows(handSelectedModels) %>% # add hand selected models
  arrange(TLG, Predator) # rearrange

# Final Models
PredictTL4_150_500 %>%
  select(-data, -knownXcol, -knownYcol, -Model, -RMSE) %>%
  rename(Pred       = Predator,
         TLGtarget = knownX) %>% 
  arrange(desc(TLG)) %>% 
  signif_df(., 3) %>% 
  kbl(align='llcccccl', caption='Final Predict. TL4 150-500 Selected Models') %>% 
  kable_paper('hover', full_width=T) %>% 
  kable_styling(fixed_thead = T)
```



## Black Bass Evaluation (Updates 2010 TMDL Staff Report Section 4.8)

### Filter for Black Bass Data
```{r}
blackBass <- fishData %>%
  filter(CommonName %in% c('Largemouth Bass', 'Smallmouth Bass', 'Spotted Bass'),
         TissueName == 'Fillet',
         Year >= 2000)
```


### QC Black Bass Data
```{r}
blackBass %>% 
  group_by(Subarea, Year) %>%
  summarise(minLength = min(TLAvgLength),
            maxLength = max(TLAvgLength),
            n = n(),
            Species = paste(unique(CommonName), collapse=', '),
            .groups='drop') %>% 
  datatable()
```


### Regression 1, Alternative 1: Standarize BB Hg Conc. to 350mm by Pooled Subareas
```{r, results="hide"}
BlkBassStd350_SubareaModels <- blackBass %>% 
  group_by(Subarea) %>%
  filter(n() > 1) %>%
  #If prediction requires extrapolation, make log model first (so selected with distinct(). Under ideal/optimized eating conditions a linear model would be appropriate but assume less than ideal environmental consumption conditions such that MeHg accumulation rate decreases as fish gets bigger/older.
  group_modify(~ predictionModels(.x, .knownYcol=Result, .knownXcol=TLAvgLength, knownX=350, preferExtrapModel='log')) %>%
  rename(Stdz350 = predictedY) %>%
  ungroup
```

#### Show Tabulated & Graphed Regression Results
<b>`r paste(nrow(BlkBassStd350_SubareaModels), 'Total Regressions Generated ~ distinct(Subarea) would select', nrow(BlkBassStd350_SubareaModels %>% distinct(Subarea)), 'regressions (1 for each Subarea)' )`</b>
```{r}
# Check Suspect Subarea GAM Regressions (identified by filter() conditions)
BlkBassStd350_SubareaModels%>%
  select(Subarea, ModelName, Stdz350, PredictionType, n) %>%
  distinct(Subarea, .keep_all=T) %>% 
  filter(grepl("gam", ModelName, ignore.case=T)) %>% 
  arrange(Subarea) %>% 
  kbl(align='llcc', caption=paste('Check These', nrow(.), 'Suspect Subarea Regressions (identified by filter() conditions)')) %>% 
  kable_paper('hover', full_width=T) %>% 
  kable_styling(fixed_thead = T)
```

<em>Model with lowest SER for each TLG & Pred Group is displayed first. Board staff manually checked that the selected models are appropriate. Board staff also checked if keeping extrapolated value is appropriate.</em><br><br>
```{r, echo=FALSE}
ui <- fluidPage(
  fluidRow(
    column(width=4, 
           selectInput('subareaName', label='Select Subarea', choices=sort(unique(BlkBassStd350_SubareaModels$Subarea), decreasing=F), selected=1)
           ),
  column(width=4,
         selectInput('subareaModel', label='Select Model for All Subarea Data', choices=NULL)
         )
  ),
  plotlyOutput('BB350Plot_subarea')
)

server <- function(input, output, session) {
  # Subarea Model Data
  data1_subarea <- reactive({
    BlkBassStd350_SubareaModels %>%
      filter(Subarea == input$subareaName)
    })
  observeEvent(data1_subarea(), {
    updateSelectInput(session, 'subareaModel',
                      choices=data1_subarea()$ModelName, selected=data1_subarea()$ModelName[1])
    })
  
  data2_subarea <- reactive({
    req(input$subareaModel)
    data1_subarea() %>%
      filter(ModelName == input$subareaModel)
    })
  
  output$BB350Plot_subarea <- renderPlotly({
    req(input$subareaModel)

    title  <- paste0('All ', input$subareaName, ' Data; ', input$subareaModel)
    xTitle <- 'Length, mm'
    yTitle <- 'Total Hg, mg/kg'
    knownXlabel <- "Std. Length"

    predictionModelPlot(data2_subarea(), knownXlabel=knownXlabel,
                        .predictedYcol=Stdz350, .groupByCol=TMDL,
                        title=title, xTitle=xTitle, yTitle=yTitle, showLegend=F)
    })
}

shinyApp(ui, server, options=list(height=500))
```

#### Select Models for Black Bass Std. 350mm Pooled Subareas
```{r}
# Subarea Models
handSelectedModels <- BlkBassStd350_SubareaModels %>% 
  filter(Subarea == 'Sacramento River'  & ModelName == 'lm' |
         Subarea == 'West Delta'        & ModelName == 'nls' |
         Subarea == 'YB-North' & ModelName == 'log')

# Replace select lowest SER models with hand selected models
BlkBassStd350_Subarea <- BlkBassStd350_SubareaModels %>%
  distinct(Subarea, .keep_all=T) %>% # select top model which should have best (i.e., smallest) SER
  anti_join(., handSelectedModels, by=c('Subarea')) %>% # remove models that will be replaced by hand selection
  bind_rows(handSelectedModels) %>% # add hand selected models
  arrange(Subarea) # rearrange 
```

#### View Pooled Subareas BB Hg Std. 350mm Hg Concentrations
```{r, fig.height=5}
# Final Models
BlkBassStd350_Subarea %>%
  select(-data, -starts_with("known"), -Model, -RMSE) %>%
  arrange(Subarea) %>% 
   mutate(Stdz350 = signif(Stdz350, 3),
         SER      = signif(SER, 3)) %>%
  kbl(align='lcccccl', caption='Final Black Bass Std 350mm Values & Models') %>% 
  kable_paper('hover', full_width=T) %>% 
  kable_styling(fixed_thead = T)
```


### Regression 1, Alternative 2: Standarize BB Hg Conc. to 350mm by Subarea & Year
```{r, results="hide"}
BlkBassStd350_YearModels <- blackBass %>% 
  group_by(Subarea, Year, Unit) %>%
  filter(n() > 1) %>%
  #If prediction requires extrapolation, make log model first (so selected with distinct(). Under ideal/optimized eating conditions a linear model would be appropriate but assume less than ideal environmental consumption conditions such that MeHg accumulation rate decreases as fish gets bigger/older.
  group_modify(~ predictionModels(.x, .knownYcol=Result, .knownXcol=TLAvgLength, knownX=350, preferExtrapModel='log')) %>%
  rename(Stdz350 = predictedY) %>%
  ungroup
```

#### Show Tabulated & Graphed Regression Results
<b>`r paste(nrow(BlkBassStd350_YearModels), 'Total Regressions Generated ~ distinct(Subarea, Year) would select', nrow(BlkBassStd350_YearModels %>% distinct(Subarea, Year)), 'regressions (1 for each Subarea & Year combination)' )`</b>
```{r}
BlkBassStd350_YearModels%>%
  select(Subarea, Year, ModelName, Stdz350, PredictionType, n) %>%
  distinct(Subarea, Year, .keep_all=T) %>% 
  filter(grepl("gam", ModelName, ignore.case=T)) %>% 
  arrange(Subarea, Year) %>% 
  kbl(align='llcc', caption=paste('Check These', nrow(.), 'Suspect Subarea by Year Regressions (identified by filter() conditions)')) %>% 
  kable_paper('hover', full_width=T) %>% 
  kable_styling(fixed_thead = T)
```
<em>Model with lowest SER for each TLG & Pred Group is displayed first. Board staff manually checked that the selected models are appropriate. Board staff also checked if keeping extrapolated value is appropriate.</em><br><br>
```{r, echo=FALSE}
ui <- fluidPage(
  fluidRow(
    column(width=4, 
           selectInput('subareaName', label='Select Subarea', choices=sort(unique(BlkBassStd350_YearModels$Subarea), decreasing=F), selected=1)
           ),
    column(width=4,
          selectInput('year', label='Select Year', choices=NULL)
          ),
    column(width=3,
         selectInput('yearModel', label='Select Model', choices=NULL)
         )
  ),
  plotlyOutput('BB350Plot_year'),
)

server <- function(input, output, session) {
  # Subarea & Year Model Data
  data1_year <- reactive({
    BlkBassStd350_YearModels %>%
      filter(Subarea == input$subareaName)
    })
  observeEvent(data1_year(), {
    updateSelectInput(session, 'year',
                      choices=data1_year()$Year, selected=data1_year()$Year[1])
    })
  
  data2_year <- reactive({
    req(input$year)
    data1_year() %>%
      filter(Year == input$year)
    })
  observeEvent(data2_year(), {
    updateSelectInput(session, 'yearModel',
                      choices=data2_year()$ModelName, selected=data2_year()$ModelName[1])
    })
  
  data3_year <- reactive({
    req(input$yearModel)
    data2_year() %>%
      filter(ModelName == input$yearModel)
    })
  
  output$BB350Plot_year <- renderPlotly({
    req(input$yearModel)
    
    title  <- paste0(input$year, ' ', input$subareaName, '; ', input$yearModel)
    xTitle <- 'Length, mm'
    yTitle <- 'Total Hg, mg/kg'
    knownXlabel <- "Std. Length"
    
    predictionModelPlot(data3_year(), knownXlabel=knownXlabel,
                        .predictedYcol=Stdz350, .groupByCol=TMDL,
                        title=title, xTitle=xTitle, yTitle=yTitle, showLegend=F)
    })
}

shinyApp(ui, server, options=list(height=500))
```

#### Select Models for Black Bass Std. 350mm by Subarea & Year
```{r}
# Subarea by Year Models
handSelectedModels <- BlkBassStd350_YearModels %>% 
  filter(Subarea == 'Central Delta' & Year == '2000' & ModelName == 'nls' |
         Subarea == 'Central Delta' & Year == '2005' & ModelName == 'gam_k=4' |
         Subarea == 'Central Delta' & Year == '2017' & ModelName == 'lm' |
         Subarea == 'Moke/Cos Rivers' & Year == '2005' & ModelName == 'nls' |
         Subarea == 'Moke/Cos Rivers' & Year == '2006' & ModelName == 'log' |
         Subarea == 'Moke/Cos Rivers' & Year == '2011' & ModelName == 'gam_k=1' |
         Subarea == 'Moke/Cos Rivers' & Year == '2017' & ModelName == 'lm' |
         Subarea == 'Moke/Cos Rivers' & Year == '2019' & ModelName == 'log' |
         Subarea == 'Sacramento River' & Year == '2005' & ModelName == 'nls' |
         Subarea == 'San Joaquin River' & Year == '2005' & ModelName == 'gam_k=1' |
         Subarea == 'San Joaquin River' & Year == '2007' & ModelName == 'log' |
         Subarea == 'West Delta' & Year == '2007' & ModelName == 'exp' |
         Subarea == 'West Delta' & Year == '2018' & ModelName == 'log' |
         Subarea == 'YB-South' & Year == '2005' & ModelName == 'log')

# Replace select lowest SER models with hand selected models
BlkBassStd350_Year <- BlkBassStd350_YearModels %>%
  distinct(Subarea, Year, .keep_all=T) %>% # select top model which should have best (i.e., smallest) SER
  anti_join(., handSelectedModels, by=c('Subarea', 'Year')) %>% # remove models that will be replaced by hand selection
  bind_rows(handSelectedModels) %>% # add hand selected models
  arrange(Subarea, Year) # rearrange 
```

#### Remove Potentially Suspect Annual Models Part I
```{r}
BlkBassStd350_YearModels%>%
  select(Subarea, Year, ModelName, Stdz350, PredictionType, n) %>%
  distinct(Subarea, Year, .keep_all=T) %>% 
  filter(n < 6 |
         PredictionType == "Extrapolated") %>% 
  arrange(Subarea, Year) %>% 
  kbl(align='llcc', caption=paste('Check These', nrow(.), 'Suspect Regressions (identified by filter() conditions)')) %>% 
  kable_paper('hover', full_width=T) %>% 
  kable_styling(fixed_thead = T)

# Sacramento River ~ 2002
# Remove: Extrapolated and Negative
BlkBassStd350_Year <- BlkBassStd350_Year %>% 
  filter(Subarea != 'Sacramento River' | Year != '2002')
```

#### Remove Potentially Suspect Annual Models Part II
```{r, fig.height=5}
BlkBassStd350_Year %>% 
  ggplot(aes(y=Stdz350, x=Subarea)) +
  geom_boxplot() +
  geom_text(aes(label=Year), nudge_x = 0.5) +
  geom_jitter(aes(colour=PredictionType), size=3) +
  theme_light() +
  theme(panel.grid.minor.y = element_blank())

# Sacramento River ~ 2001
# Remove: Extreme Extrapolation
BlkBassStd350_Year <- BlkBassStd350_Year %>% 
  filter(Subarea != 'Sacramento River' | Year != '2001')

# San Joaquin River ~ 2018
# Remove: Extreme Extrapolation
BlkBassStd350_Year <- BlkBassStd350_Year %>% 
  filter(Subarea != 'San Joaquin River' | Year != '2018')

# Yolo Bypass North ~ 2003
# Remove: Extreme Extrapolation
BlkBassStd350_Year <- BlkBassStd350_Year %>% 
  filter(Subarea != 'YB-North' | Year != '2003')
```

#### View Subarea by Year BB Hg Std. 350mm Hg Concentrations
```{r, fig.height=5}
BlkBassStd350_Year %>% 
  ggplot(aes(y=Stdz350, x=Subarea)) +
  geom_boxplot() +
  geom_text(aes(label=Year), nudge_x = 0.5) +
  geom_jitter(aes(colour=PredictionType), size=3) +
  theme_light() +
  theme(panel.grid.minor.y = element_blank())

# Final Models
BlkBassStd350_Year %>%
  select(-data, -starts_with("known"), -Model, -RMSE) %>%
  arrange(Subarea, Year) %>% 
   mutate(Stdz350 = signif(Stdz350, 3),
         SER      = signif(SER, 3)) %>%
  kbl(align='lcccccl', caption='Final Black Bass Std 350mm Values & Models') %>% 
  kable_paper('hover', full_width=T) %>% 
  kable_styling(fixed_thead = T)
```


### Compare Alternative 1 and 2: Black Bass Std. 350mm by Subarea & Year vs. Pooled Subareas
```{r}
# Because NumberFishPerComp is not carried over in BlkBassStd350_Year, create df that includes sums(NumberFishPerComp) in same grouping
totalFish_SubareaYear <- blackBass %>% 
  mutate(Year = year(SampleDate)) %>% 
  group_by(Subarea, Year) %>% 
  summarise(Std350Samples = n(),                      # create Std350Samples column
            Std350TotalFish = sum(NumberFishPerComp), # create Std350TotalFish column
            .groups = 'drop')

# Calc yearly wieghted avg for each subarea
BlkBassStd350_Subarea.weightedYear <- BlkBassStd350_Year %>% 
  left_join(., totalFish_SubareaYear, by=c('Subarea', 'Year')) %>%  # add Std350Samples & Std350TotalFish column
  group_by(Subarea) %>%
  summarise(Stdz350       = weighted.mean(Stdz350, Std350TotalFish),
            Std350Samples = sum(Std350Samples),
            Std350TotalFish = sum(Std350TotalFish),
            .groups = 'drop') %>% 
  mutate(Method = 'Yearly Wt. Avg.')


# Each Subarea Pooled Std 350 Hg Regression
BlkBassStd350_Subarea.pooled <- BlkBassStd350_Subarea %>%
  select(Subarea, Stdz350, n) %>% 
  rename(Std350TotalFish = n) %>% 
  mutate(Std350Samples = NA,
         Method = 'Pooled Regression')

# Compare Results
BlkBassStd350_Subarea.avgHg <- bind_rows(BlkBassStd350_Subarea.weightedYear, BlkBassStd350_Subarea.pooled) %>% 
  # mutate(`_` = 'Black Bass Hg Std. 350mm') %>%
  select(Subarea, Method, Stdz350) %>%
  pivot_wider(names_from=c(Method), values_from=c(Stdz350))

BlkBassStd350_Subarea.avgHg %>% 
  round_df(., 3) %>% 
  kbl(caption='Reproduce TMDL Staff Report Table 4.10') %>% 
  kable_paper('hover', full_width=F)

writexl::write_xlsx(BlkBassStd350_Subarea.avgHg,
                    paste0(wd, '/Reeval_Impl_Goals_Linkage_Analysis/eval1_TLG Eval & BB Impl Goal/Output/TLG Appdx Table X.X 2000-2019 BB Std 350 Annual WtAvg vs Pooled.xlsx'))
```

### Regression 2, Alternative 1: BB Implementation Goal by Pooled Subareas

#### Add Annual Black Bass 350mm Stand. Length Avg. Results to TLG Hg Composite Wt. Avg Results
```{r}
TLG_BlkBass_Subarea.pooled <- TLG_AvgHg_subarea %>% #TLG_AvgHg_subarea %>% 
  left_join(., BlkBassStd350_Subarea.pooled %>% select(Subarea, Stdz350, Std350Samples, Std350TotalFish), by='Subarea') %>%  # adds Stdz350 & Std350Samples columns
  filter(!is.na(Stdz350))
```

#### Perform Regressions to Determine Black Bass Std. 350mm Predicted Safe Concentrations
```{r, results="hide"}
BBSafeLevels_SubareaPoolModels <- TLG_BlkBass_Subarea.pooled %>%  
  group_by(TLG, Predator) %>%
  group_modify(~ predictionModels(.x, .knownYcol=Stdz350, .knownXcol=CompositeWeightedAvg, knownX=.x$TLGtarget[1], preferExtrapModel=NULL)) %>%
  rename(PredictBBSafeLevel = predictedY) %>% 
  ungroup
```

#### View Black Bass Std. 350mm Predicted Safe Concentration Regressions
<b>`r paste(nrow(BBSafeLevels_SubareaPoolModels), 'Total Regressions Generated ~ distinct(TLG) would select', nrow(BBSafeLevels_SubareaPoolModels %>% distinct(TLG)), 'regressions (1 for each TLG)' )`</b>
```{r}
# Check Suspect Regressions (identified by filter() conditions)
BBSafeLevels_SubareaPoolModels%>%
  select(TLG, Predator, ModelName, PredictBBSafeLevel, PredictionType, n) %>%
  distinct(TLG, PredictionType, .keep_all=T) %>% 
  filter(grepl("gam", ModelName, ignore.case=T) |
           PredictionType == "Extrapolated") %>%
  arrange(desc(TLG)) %>% 
  kbl(align='llcc', caption=paste('Check These', nrow(.), 'Suspect Regressions (identified by filter() conditions)')) %>% 
  kable_paper('hover', full_width=T) %>% 
  kable_styling(fixed_thead = T)
```
<em>Model with lowest SER for each TLG & Pred Group is displayed first. gam models tend to have lowest SER, but need to check that model makes an appropriate regression. Also need to check if keeping extrapolated value is appropriate.</em><br><br>
```{r, echo=FALSE}
ui <- fluidPage(
  fluidRow(
    column(width=3, 
           selectInput('tlgName', label='Select X-axis TLG', choices=sort(unique(BBSafeLevels_SubareaPoolModels$TLG), decreasing=T), selected=1)
           ),
  column(width=3,
         selectInput('predName', label='Select Predatory Group', choices=NULL)
         ),
  column(width=3,
         selectInput('model', label='Select Model', choices=NULL)
         )
  ),
  plotlyOutput('tlgPlot')
)

server <- function(input, output, session) {
  data1 <- reactive({
    BBSafeLevels_SubareaPoolModels %>%
      filter(TLG == input$tlgName)
    })
  observeEvent(data1(), {
    updateSelectInput(session, 'predName',
                      choices=data1()$Predator, selected=data1()$Predator[1])
    })
  
  data2 <- reactive({
    req(input$predName)
    data1() %>%
      filter(Predator == input$predName)
    })
   observeEvent(data2(), {
    updateSelectInput(session, 'model',
                      choices=data2()$ModelName, selected=data2()$ModelName[1])
    })
   
   data3 <- reactive({
    req(input$model)
    data2() %>%
      filter(ModelName == input$model)
   })
  
  output$tlgPlot <- renderPlotly({
    req(input$model)
    
    title  <- paste('Std.350 BB vs.', data3()$TLG[1])
    xTitle <- paste(data3()$TLG[1],'(mg/kg)')
    yTitle <- 'Std.350 BB (mg/kg)'
    knownXlabel <- paste0("TLG Target<sub>",data2()$Predator,"</sub>")
    
    predictionModelPlot(data3(), knownXlabel=knownXlabel,
                    .predictedYcol=PredictBBSafeLevel, .groupByCol=Subarea, allGroups=sort(unique(fishData$Subarea)),
                    title=title, xTitle=xTitle, yTitle=yTitle, showLegend=T)
    })
}

shinyApp(ui, server, options=list(height=500))
```

#### Select Models for Black Bass Std. 350mm Predicted Safe Concentration
```{r}
handSelectedModels <- BBSafeLevels_SubareaPoolModels %>% 
  filter(TLG == 'TL3 Fish (50 - <150 mm)' & ModelName == 'log')

# Replace select lowest SER models with hand selected models
BBSafeLevels_SubareaPool <- BBSafeLevels_SubareaPoolModels %>%
  distinct(TLG, Predator, .keep_all=T) %>% # select top model which should have best (i.e., smallest) SER
  anti_join(., handSelectedModels, by=c('TLG')) %>% # remove models that will be replaced by hand selection
  bind_rows(handSelectedModels) %>% # add hand selected models
  arrange(desc(TLG), Predator) # rearrange 

# Final Selected Models
BBSafeLevels_SubareaPool%>%
  select(-data, -starts_with("known"), -Model) %>%
  signif_df(., 3) %>% 
  kbl(caption='Final Black Bass Safe Level Selected Models') %>% 
  kable_paper('hover', full_width=F)

# Avg SER of all Models
avgSER_SubareaPool <- BBSafeLevels_SubareaPool %>%
  distinct(TLG, SER) %>% 
  summarise(mean(SER)) %>% pull()
avgSER_SubareaPool
```

#### Reproduce Table 4.9
Reproduce 2010 TMDL Staff Report Table 4.9 (Predicted Safe Concentrations of Methylmercury in 150-500 mm TL4 Fish
and Standard 350-mm Largemouth Bass Corresponding to Trophic Level
Food Group (TLFG) Targets for the Protection of Piscivorous Species)
```{r, echo=FALSE}
TL4targets <- PredictTL4_150_500 %>% 
  select(TLG, Predator, PredictTL4_150_500)

BBSafeLevels_SubareaPool %>% 
  select(TLG, Predator, knownX, PredictBBSafeLevel) %>% 
  left_join(., TL4targets, by=c('TLG', 'Predator')) %>% 
  rename(Target = knownX) %>% 
  select(TLG, Predator, Target, PredictTL4_150_500, PredictBBSafeLevel) %>% 
  arrange(desc(TLG), Target) %>% 
  signif_df(., 3) %>% 
  mutate(PredictTL4_150_500 = case_when(is.na(PredictTL4_150_500) ~ Target,
                                      T ~ PredictTL4_150_500)) %>% 
  kbl(format="html", escape=T) %>% 
  kable_paper('hover', full_width=F)
```

#### Black Bass Standard 350mm Hg Implementation Goal
```{r}
# Select Lowest PredictBBSafeLevel as Imp. Goal
BBSafeLevel_Pool_final <- BBSafeLevels_SubareaPool %>% 
  group_by(knownYcol) %>% 
  filter(PredictBBSafeLevel == min(PredictBBSafeLevel))

title  <- 'Black Bass Imp. Goal using Pooled Black Bass Std. 350mm [Hg]'
xTitle <- paste(BBSafeLevel_Pool_final$TLG[1],'(mg/kg)')
yTitle <- 'Std.350 BB (mg/kg)'
ImplGoalXlabel <- paste0("TLG Target<sub>",BBSafeLevel_Pool_final$Predator,"</sub>")

predictionModelPlot(BBSafeLevel_Pool_final, knownXlabel=ImplGoalXlabel,
                    .predictedYcol=PredictBBSafeLevel, .groupByCol=Subarea, allGroups=sort(unique(fishData$Subarea)),
                    title=title, xTitle=xTitle, yTitle=yTitle, showLegend=T)
```


### Regression 2, Alternative 2: BB Implementation Goal Weighted by Sample Size for Each Subarea
Calculate Annual Black Bass 350mm Stand. Length Avg. Weighted by Sample Size for Each Subarea

#### Add Annual Black Bass 350mm Stand. Length Avg. Results to TLG Hg Composite Wt. Avg Results
```{r}
TLG_BlkBass_Subarea.weightedYear <- TLG_AvgHg_subarea %>% #TLG_AvgHg_subarea %>% 
  left_join(., BlkBassStd350_Subarea.weightedYear %>% select(Subarea, Stdz350, Std350Samples, Std350TotalFish), by='Subarea') %>%  # adds Stdz350 & Std350Samples columns
  filter(!is.na(Stdz350))
```

#### Perform Regressions to Determine Black Bass Std. 350mm Predicted Safe Concentrations
```{r, results="hide"}
BBSafeLevels_SubareaWtYrModels <- TLG_BlkBass_Subarea.weightedYear %>%  
  group_by(TLG, Predator) %>%
  group_modify(~ predictionModels(.x, .knownYcol=Stdz350, .knownXcol=CompositeWeightedAvg, knownX=.x$TLGtarget[1], preferExtrapModel=NULL)) %>%
  rename(PredictBBSafeLevel = predictedY) %>% 
  ungroup
```

#### View Black Bass Std. 350mm Predicted Safe Concentration Regressions
<b>`r paste(nrow(BBSafeLevels_SubareaWtYrModels), 'Total Regressions Generated ~ distinct(TLG) would select', nrow(BBSafeLevels_SubareaWtYrModels %>% distinct(TLG)), 'regressions (1 for each TLG)' )`</b>
```{r}
# Check Suspect Regressions (identified by filter() conditions)
BBSafeLevels_SubareaWtYrModels%>%
  select(TLG, Predator, ModelName, PredictBBSafeLevel, PredictionType, n) %>%
  distinct(TLG, PredictionType, .keep_all=T) %>% 
  filter(grepl("gam", ModelName, ignore.case=T) |
           PredictionType == "Extrapolated") %>%
  arrange(desc(TLG)) %>% 
  kbl(align='llcc', caption=paste('Check These', nrow(.), 'Suspect Regressions (identified by filter() conditions)')) %>% 
  kable_paper('hover', full_width=T) %>% 
  kable_styling(fixed_thead = T)
```
<em>Model with lowest SER for each TLG & Pred Group is displayed first. gam models tend to have lowest SER, but need to check that model makes an appropriate regression. Also need to check if keeping extrapolated value is appropriate.</em><br><br>
```{r, echo=FALSE}
ui <- fluidPage(
  fluidRow(
    column(width=3, 
           selectInput('tlgName', label='Select X-axis TLG', choices=sort(unique(BBSafeLevels_SubareaWtYrModels$TLG), decreasing=T), selected=1)
           ),
  column(width=3,
         selectInput('predName', label='Select Predatory Group', choices=NULL)
         ),
  column(width=3,
         selectInput('model', label='Select Model', choices=NULL)
         )
  ),
  plotlyOutput('tlgPlot')
)

server <- function(input, output, session) {
  data1 <- reactive({
    BBSafeLevels_SubareaWtYrModels %>%
      filter(TLG == input$tlgName)
    })
  observeEvent(data1(), {
    updateSelectInput(session, 'predName',
                      choices=data1()$Predator, selected=data1()$Predator[1])
    })
  
  data2 <- reactive({
    req(input$predName)
    data1() %>%
      filter(Predator == input$predName)
    })
   observeEvent(data2(), {
    updateSelectInput(session, 'model',
                      choices=data2()$ModelName, selected=data2()$ModelName[1])
    })
   
   data3 <- reactive({
    req(input$model)
    data2() %>%
      filter(ModelName == input$model)
   })
  
  output$tlgPlot <- renderPlotly({
    req(input$model)
    
    title  <- paste('Std.350 BB vs.', data3()$TLG[1])
    xTitle <- paste(data3()$TLG[1],'(mg/kg)')
    yTitle <- 'Std.350 BB (mg/kg)'
    knownXlabel <- paste0("TLG Target<sub>",data2()$Predator,"</sub>")
    
    predictionModelPlot(data3(), knownXlabel=knownXlabel,
                    .predictedYcol=PredictBBSafeLevel, .groupByCol=Subarea, allGroups=sort(unique(fishData$Subarea)),
                    title=title, xTitle=xTitle, yTitle=yTitle, showLegend=T)
    })
}

shinyApp(ui, server, options=list(height=500))
```

#### Select Models for Black Bass Std. 350mm Predicted Safe Concentration
```{r}
handSelectedModels <- BBSafeLevels_SubareaWtYrModels %>% 
  filter(TLG == 'none' & ModelName == 'none')

# Replace select lowest SER models with hand selected models
BBSafeLevels_SubareaWtYr <- BBSafeLevels_SubareaWtYrModels %>%
  distinct(TLG, Predator, .keep_all=T) %>% # select top model which should have best (i.e., smallest) SER
  anti_join(., handSelectedModels, by=c('TLG')) %>% # remove models that will be replaced by hand selection
  bind_rows(handSelectedModels) %>% # add hand selected models
  arrange(desc(TLG), Predator) # rearrange 

# Final Selected Models
BBSafeLevels_SubareaWtYr %>%
  select(-data, -starts_with("known"), -Model) %>%
  signif_df(., 3) %>% 
  kbl(caption='Final Black Bass Safe Level Selected Models') %>% 
  kable_paper('hover', full_width=F)

# Avg SER of all Models
avgSER_SubareaWtYr <- BBSafeLevels_SubareaWtYr %>%
  distinct(TLG, SER) %>% 
  summarise(mean(SER)) %>% pull()
avgSER_SubareaWtYr

```

#### Reproduce Table 4.9
Reproduce 2010 TMDL Staff Report Table 4.9 (Predicted Safe Concentrations of Methylmercury in 150-500 mm TL4 Fish
and Standard 350-mm Largemouth Bass Corresponding to Trophic Level
Food Group (TLFG) Targets for the Protection of Piscivorous Species)
```{r, echo=FALSE}
TL4targets <- PredictTL4_150_500 %>% 
  select(TLG, Predator, PredictTL4_150_500)
Table_4.9 <- BBSafeLevels_SubareaWtYr %>% 
  select(TLG, Predator, knownX, PredictBBSafeLevel) %>% 
  left_join(., TL4targets, by=c('TLG', 'Predator')) %>% 
  rename(Target = knownX) %>% 
  select(TLG, Predator, Target, PredictTL4_150_500, PredictBBSafeLevel) %>% 
  arrange(desc(TLG), Target) %>% 
  signif_df(., 3) %>% 
  mutate(PredictTL4_150_500 = case_when(is.na(PredictTL4_150_500) ~ Target,
                                      T ~ PredictTL4_150_500)) %>% 
  rename(`Trophic Level Group` = TLG)

Table_4.9 %>% 
  kbl(format="html", escape=T) %>% 
  kable_paper('hover', full_width=F)
```

#### Black Bass Standard 350mm Hg Implementation Goal
```{r}
# Select Lowest PredictBBSafeLevel as Imp. Goal
BBSafeLevel_WtYr_final <- BBSafeLevels_SubareaWtYr %>% 
  group_by(knownYcol) %>% 
  filter(PredictBBSafeLevel == min(PredictBBSafeLevel))

title  <- 'Black Bass Imp. Goal using Wt Year Avg. Black Bass Std. 350mm [Hg]'
xTitle <- paste(BBSafeLevel_WtYr_final$TLG[1],'(mg/kg)')
yTitle <- 'Std.350 BB (mg/kg)'
ImplGoalXlabel <- paste0("TLG Target<sub>",BBSafeLevel_WtYr_final$Predator,"</sub>")

predictionModelPlot(BBSafeLevel_WtYr_final, knownXlabel=ImplGoalXlabel,
                    .predictedYcol=PredictBBSafeLevel, .groupByCol=Subarea, allGroups=sort(unique(fishData$Subarea)),
                    title=title, xTitle=xTitle, yTitle=yTitle, showLegend=T)
```




## Final Black Bass Implementation Goal for Linkage Analysis
Selected Yearly Weighted: Black Bass Std. 350mm Hg Imp. Goal of <b>`r signif(BBSafeLevel_WtYr_final %>% pull(PredictBBSafeLevel), 3)`</b> as the final value to be used in the linkage analysis. The Black Bass Implementation Goal regressions using yearly weighted Black Bass 350mm Std Hg concentrations had a lower average SER of `r signif(avgSER_SubareaWtYr, 3)` mg/kg vs. `r signif(avgSER_SubareaPool, 3)` mg/kg for the pooled Black Bass 350mm Std Hg concentrations. Also, the North Yolo Bypass did not have sufficient black bass data (n=3) to be as representative as other subareas, so Board staff decided to exclude this subarea from the Black Bass Implementation Goal regressions.

## Save Data used for Evaluations
```{r}
fishData_save <- fishData %>% 
  select(Subarea, SourceID, ProjectCode, StationName, SampleDate, Analyte, Unit, Result, MDL, RL, ResultQualCode, Censored, CommonName, TLAvgLength, TissueName, NumberFishPerComp) %>%
  rename(ReferenceCode=ProjectCode,
         TLAvgLength_mm = TLAvgLength) %>% 
  mutate(SourceID = case_when(SourceID == "CEDENFish" ~ "CEDEN",
                              T ~ SourceID),
         ReferenceCode = case_when(SourceID == "R5MF" ~ ReferenceCode,
                                   T ~ NA_character_)) %>% 
  arrange(Subarea, SampleDate, StationName)

blackBass_save <- blackBass %>% 
  select(Subarea, SourceID, ProjectCode, StationName, SampleDate, Analyte, Unit, Result, MDL, RL, ResultQualCode, Censored, CommonName, TLAvgLength, TissueName, NumberFishPerComp) %>%
  rename(ReferenceCode=ProjectCode,
         TLAvgLength_mm = TLAvgLength) %>% 
  mutate(SourceID = case_when(SourceID == "CEDENFish" ~ "CEDEN",
                              T ~ SourceID),
         ReferenceCode = case_when(SourceID == "R5MF" ~ ReferenceCode,
                                   T ~ NA_character_)) %>% 
  arrange(Subarea, SampleDate, StationName)

writexl::write_xlsx(list("Data Considered for TLG Eval"=fishData_save,
                         "Data Considered for BB Eval"=blackBass_save),
                    paste0(wd, '/Reeval_Impl_Goals_Linkage_Analysis/eval1_TLG Eval & BB Impl Goal/Output/Appdx A_Data Compilation_Data Considered for TLG & BB Eval.xlsx'))
```
## Error: <text>:2:8: unexpected '|'
## 1: ---
## 2: title: |
##           ^

The R session information (including the OS info, R version and all packages used):

    sessionInfo()
## R version 4.2.2 (2022-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 22621)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_United States.utf8  LC_CTYPE=English_United States.utf8   
## [3] LC_MONETARY=English_United States.utf8 LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.utf8    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] lubridate_1.8.0    plotly_4.10.0      readxl_1.4.1       actuar_3.3-0      
##  [5] NADA_1.6-1.1       forcats_0.5.2      stringr_1.4.1      dplyr_1.0.9       
##  [9] purrr_0.3.4        readr_2.1.2        tidyr_1.2.0        tibble_3.1.8      
## [13] ggplot2_3.3.6      tidyverse_1.3.2    fitdistrplus_1.1-8 survival_3.4-0    
## [17] MASS_7.3-58.1     
## 
## loaded via a namespace (and not attached):
##  [1] lattice_0.20-45     assertthat_0.2.1    digest_0.6.29       utf8_1.2.2         
##  [5] R6_2.5.1            cellranger_1.1.0    backports_1.4.1     reprex_2.0.2       
##  [9] evaluate_0.16       highr_0.9           httr_1.4.4          pillar_1.8.1       
## [13] rlang_1.0.5         lazyeval_0.2.2      googlesheets4_1.0.1 rstudioapi_0.14    
## [17] data.table_1.14.2   Matrix_1.5-1        splines_4.2.2       googledrive_2.0.0  
## [21] htmlwidgets_1.5.4   munsell_0.5.0       broom_1.0.1         compiler_4.2.2     
## [25] modelr_0.1.9        xfun_0.32           pkgconfig_2.0.3     htmltools_0.5.3    
## [29] tidyselect_1.1.2    fansi_1.0.3         viridisLite_0.4.1   crayon_1.5.1       
## [33] tzdb_0.3.0          dbplyr_2.2.1        withr_2.5.0         grid_4.2.2         
## [37] jsonlite_1.8.0      gtable_0.3.1        lifecycle_1.0.1     DBI_1.1.3          
## [41] magrittr_2.0.3      scales_1.2.1        writexl_1.4.0       cli_3.3.0          
## [45] stringi_1.7.8       fs_1.5.2            xml2_1.3.3          ellipsis_0.3.2     
## [49] generics_0.1.3      vctrs_0.4.1         expint_0.1-7        tools_4.2.2        
## [53] glue_1.6.2          hms_1.1.2           fastmap_1.1.0       colorspace_2.0-3   
## [57] gargle_1.2.0        rvest_1.0.3         knitr_1.40          haven_2.5.1
    Sys.time()
## [1] "2024-01-05 12:02:20 PST"