This report is automatically generated with the R
package knitr
(version 1.40
)
.
--- title: | | TLG Eval. & Black Bass Hg Imp. Goal (2002-2019 data) | Final Eval Results Using Data Pooled by Subarea & Year with Overall Mean Weighted by Yearly Sample Size 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 new data for fish species and black bass from 2002 through 2019. It shows the final results of determining the BB Implementation Goal based on the two evaluations. The rmarkdown "1.TEST TLG 1998-2019 & BB Impl Goal 2000-2019_pooled vs weighted avg.Rmd" showed that pooling data by year and subarea provided the lowest avg SER of predicted BB protective concentration models. The next evaluation for year ranges in rmarkdown "2.TEST TLG & BB_ year ranges.Rmd" showed that data years 2002-2019 provided the lowest avg SER of predicted BB protective concentration models. ```{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") ``` ## 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)) # Supporting File - Trophic Level Group (TLG) reference # This file was developed to be consistent with the 2010 TMDL Staff Report and Statewide Mercury Provisions (e.g. a juvenile Black Bass may be TL3, but an adult may be TL4) 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) # Convert columns that are character but should be numeric or date TLgroups <- chara_to_NumDate(TLgroups) ``` ## Initial Fish Data Prep ```{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 Trophic Level Group (TLG) ```{r} fishData3 <- fishData2 %>% mutate(commonname = tolower(CommonName)) %>% # create column of Trophic Level (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() # For updated revaluation, the years 2002 through 2019 (new fish data) were determined to provide the lowest average SER for Black Bass regressions (see "2.TEST TLG & BB_ year ranges.Rmd") fishData4 <- fishData3 %>% filter(Year >= 2002) ``` ### 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 fishData3 %>% 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 %>% 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 custom function 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 food groups."</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='Compare Results of Using Both Tissue Types vs. Either Fillet or Whole Body') %>% kable_paper('hover', full_width=F) %>% kable_styling(fixed_thead = T) # The table above compares 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) # Board staff evaluated fish sampling locations using GIS. No fish were sampled from the Marsh Creek subarea between 1998 to 2019. Fish sampling in the North and South Yolo Bypass between 2002 to 2019 provided more representative data for these subareas, which were included in trophic level group evaluation. ``` ## 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$gamma$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) ``` ### Finalize QC'd Data ```{r} fishData <- fishData5 ``` ## TLG Evaluation (Updates 2010 TMDL Staff Report Section 4.7) ### Calcualte TLG Composite Weighted Average for each Subarea ```{r} # Add "TL" to numeric 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)) # Total Number of Composites, Fish & Species totalFishSummary <- TLG_AvgHg_subarea %>% summarize(compositeTotal = sum(CompositeSamples), fishTotal = sum(CompositeTotalFish)) names(totalFishSummary$compositeTotal) <- "Total Number of Composite Samples" names(totalFishSummary$fishTotal) <- "Total Number of Fish" totalSpecies <- length(unique(TLG$CommonName)) names(totalSpecies) <- "Total Number of Species" totalFishSummary$compositeTotal; totalFishSummary$fishTotal; totalSpecies ``` #### 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) 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) setwd(wd) writexl::write_xlsx(Table_4.7, 'Reeval_Impl_Goals_Linkage_Analysis/eval1_TLG Eval & BB Impl Goal/Output/Table 4.7.xlsx') ``` ### Join TL4 150-500mm Results to Other TLG Results ```{r} # Prepares data for regression # 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 (150-500 mm)' & ModelName == 'log'| TLG == 'TL3 Fish (150-350 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) # Because human consumption of TL4 and TL3 fish remain the lowest predicted TL4 (150-500 mm) protective concentration of TL3 and TL4 TLGs, Board staff conclude that the WQOs proposed in the 2010 Staff Report still ensure safe mercury concentrations in Delta fish for consumption by humans and wildlife. ``` ### Save TL4 150-500 Model Datafram to create Figs of Regression Models ```{r} save(PredictTL4_150_500, file=paste0(wd, '/Reeval_Impl_Goals_Linkage_Analysis/eval1_TLG Eval & BB Impl Goal/Output/Regression Models_Predicted TL4 150-500 Protectiv Conc.RData')) ``` ## Black Bass Evaluation (Updates 2010 TMDL Staff Report Section 4.8) ### Filter for Black Bass Data ```{r} blackBass_2002 <- fishData %>% filter(CommonName %in% c('Largemouth Bass', 'Smallmouth Bass', 'Spotted Bass'), TissueName == 'Fillet', Year >= 2002) ``` ### Black Bass Data QC ```{r} blackBass_2002_dataSummary <- blackBass_2002 %>% group_by(Subarea, Year) %>% summarise(LengthRange = paste(min(TLAvgLength),"-",max(TLAvgLength)), n = n(), Species = paste(unique(CommonName), collapse=', '), .groups='drop') blackBass_2002_dataSummary %>% datatable() ``` ### Regression 1: Standarize BB Hg Conc. to 350mm by Subarea & Year ```{r, results="hide"} BlkBassStd350_Models_2002 <- blackBass_2002 %>% 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 %>% mutate(BlackBassEval = "Included") %>% # assume all std 250 values will be included in Black Bass Eval and then exclude as needed based on further evaluation below left_join(., blackBass_2002_dataSummary %>% select(Subarea, Year, LengthRange), # adds LengthRange Column by=c("Subarea", "Year")) ``` #### Show Tabulated & Graphed Regression Results <b>`r paste(nrow(BlkBassStd350_Models_2002), 'Total Regressions Generated ~ distinct(Subarea, Year) would select', nrow(BlkBassStd350_Models_2002 %>% distinct(Subarea, Year)), 'regressions (1 for each Subarea & Year combination)' )`</b> ```{r} BlkBassStd350_Models_2002 %>% select(Subarea, Year, ModelName, Stdz350, PredictionType, n, BlackBassEval) %>% 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_Models_2002$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_Models_2002 %>% 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_Models_2002 %>% filter(Subarea == 'Central Delta' & Year == '2005' & ModelName == 'gam_k=4' | Subarea == 'Central Delta' & Year == '2016' & ModelName == 'gam_k=1' | 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 == '2007' & ModelName == 'log' | Subarea == 'West Delta' & Year == '2007' & ModelName == 'exp' | Subarea == 'West Delta' & Year == '2018' & ModelName == 'log' | Subarea == 'none' & Year == 'none' & ModelName == 'none') # Replace select lowest SER models with hand selected models BlkBassStd350_SelectedModels_2002 <- BlkBassStd350_Models_2002 %>% 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_Models_2002%>% distinct(Subarea, Year, .keep_all=T) %>% filter(n < 6 | PredictionType == "Extrapolated") %>% select(Subarea, Year, ModelName, Stdz350, PredictionType, ExtrapLength, n) %>% 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 # Exclude: Extrapolated and Negative BlkBassStd350_SelectedModels_2002 <- BlkBassStd350_SelectedModels_2002 %>% mutate(BlackBassEval = case_when(Subarea == 'Sacramento River' & Year == '2002' ~ "Excluded", T ~ BlackBassEval)) ``` #### Remove Potentially Suspect Annual Models Part II ```{r, fig.height=5} BlkBassStd350_SelectedModels_boxplot <- BlkBassStd350_SelectedModels_2002 %>% filter(BlackBassEval == "Included") %>% mutate(Year=ifelse(PredictionType == "Extrapolated", Year, ""), Subarea = gsub(" River", "\nRiver", Subarea)) # let Subareas with "River" split across 2 lines colors_graphs <- RColorBrewer::brewer.pal(n = 11, name = 'RdYlBu') # "#A50026" "#D73027" "#F46D43" "#FDAE61" "#FEE090" "#FFFFBF" "#E0F3F8" "#ABD9E9" "#74ADD1" "#4575B4" "#313695" used_colors <- c("Interpolated" = colors_graphs[9], "Extrapolated"=colors_graphs[3]) BlkBassStd350_SelectedModels_boxplot %>% ggplot(aes(y=Stdz350, x=Subarea)) + geom_boxplot(data=. %>% filter(PredictionType == "Interpolated")) + geom_point(data=BlkBassStd350_SelectedModels_boxplot, aes(fill=PredictionType), size=3, colour="black",pch=21) + # pch adds black circle around data point scale_fill_manual(values = used_colors[unique(BlkBassStd350_SelectedModels_boxplot$PredictionType)]) + geom_text(aes(label=Year), nudge_x = -0.235, nudge_y = -0.018) + xlab(element_blank()) + ylab("Standardized 350 mm Black Bass Hg Concentration (mg/Kg)") + labs(fill = "Prediction Type") + theme_light() + theme(text = element_text(size=14, family="sans"), axis.title=element_text(size=12, face="bold"), axis.text.x = element_text(size=10.5, angle = 0, hjust = 0.5), panel.grid.major.y=element_line(size=1), panel.grid.minor=element_blank(), legend.text = element_text(size = 14), legend.title = element_text(size = 14, face="bold")) ggsave(filename = paste0(wd, '/Reeval_Impl_Goals_Linkage_Analysis/eval1_TLG Eval & BB Impl Goal/Output/Fig 4.2 Std350 Black Bass Concentration.png'), # plot = , width = 9, height = 5.75, units = "in", dpi = 125) # San Joaquin River ~ 2018 # Exclude: Extreme Extrapolation BlkBassStd350_SelectedModels_2002 <- BlkBassStd350_SelectedModels_2002 %>% mutate(BlackBassEval = case_when(Subarea == 'San Joaquin River' & Year == '2018' ~ "Excluded", T ~ BlackBassEval)) # Yolo Bypass North ~ 2003 # Exclude: Extreme Extrapolation BlkBassStd350_SelectedModels_2002 <- BlkBassStd350_SelectedModels_2002 %>% mutate(BlackBassEval = case_when(Subarea == 'YB-North' & Year == '2003' ~ "Excluded", T ~ BlackBassEval)) # Save Model Evaluation Results BlkBassStd350_SelectedModels_2002_save <- BlkBassStd350_SelectedModels_2002 %>% select(BlackBassEval, Subarea, everything()) %>% arrange(BlackBassEval, Subarea) writexl::write_xlsx(BlkBassStd350_SelectedModels_2002_save, paste0(wd, '/Reeval_Impl_Goals_Linkage_Analysis/eval1_TLG Eval & BB Impl Goal/Output/Table 4.4 Black Bass Model Evaluation Results.xlsx')) ``` #### View Subarea by Year BB Hg Std. 350mm Hg Concentrations ```{r, fig.height=5} BlkBassStd350_FinalModels_2002 <- BlkBassStd350_SelectedModels_2002 %>% filter(BlackBassEval == "Included") BlkBassStd350_FinalModels_boxplot <- BlkBassStd350_FinalModels_2002 %>% mutate(Year=ifelse(PredictionType == "Extrapolated", Year, ""), Subarea = gsub(" River", "\nRiver", Subarea)) # let Subareas with "River" split across 2 lines) BlkBassStd350_FinalModels_boxplot %>% ggplot(aes(y=Stdz350, x=Subarea)) + geom_boxplot(data=. %>% filter(PredictionType == "Interpolated")) + geom_point(aes(fill=PredictionType), size=3, colour="black",pch=21) + # pch adds black circle around data point geom_text(aes(label=Year), nudge_x = -0.235, nudge_y = -0.018) + xlab("") + ylab("Standardized 350 mm Black Bass Hg Concentration (mg/Kg)") + theme_light() + theme(text = element_text(size=14, family="sans"), axis.title=element_text(size=12, face="bold"), axis.text.x = element_text(size=10.5, angle = 0, hjust = 0.5), panel.grid.major.y=element_line(size=1), panel.grid.minor=element_blank(), legend.text = element_text(size = 14), legend.title = element_text(size = 14, face="bold")) # Final Models BlkBassStd350_FinalModels_2002 %>% 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) # Combine All Fish & Final BB Sampling Locations BBdata <- NULL for(i in 1:nrow(BlkBassStd350_FinalModels_2002)){ BBdata <- BBdata %>% bind_rows(BlkBassStd350_FinalModels_2002$data[[i]] %>% mutate(Subarea = BlkBassStd350_FinalModels_2002$Subarea[[i]])) } bbSampLoc <- BBdata %>% mutate(Year = year(SampleDate), Species = 'Black Bass', 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) fish_bb_SampLoc <- fishData4 %>% filter(!grepl("largemouth|smallmouth|spotted", CommonName, ignore.case=T)) %>% # remove black bass species so these locations are not shown as "All Species" when shown by Black Bass alone on sampling location map 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) %>% bind_rows(bbSampLoc) setwd(wd) writexl::write_xlsx(fish_bb_SampLoc, path='Reeval_Impl_Goals_Linkage_Analysis/eval1_TLG Eval & BB Impl Goal/Output/TLG Eval_Fish&BB_SampleLocations.xlsx') ``` #### Calculate Annual Black Bass 350mm Stand. Length Avg. Weighted by Sample Size for Each Subarea ```{r} # Because NumberFishPerComp is not carried over in BlkBassStd350_FinalModels_2002, create data frame that includes sums(NumberFishPerComp) in same grouping totalFish_SubareaYear <- blackBass_2002 %>% mutate(Year = year(SampleDate)) %>% group_by(Subarea, Year) %>% summarise(Std350Samples = n(), # create Std350Samples column Std350TotalFish = sum(NumberFishPerComp), # create Std350TotalFish column .groups = 'drop') # Calculate annual weighted avg for each subarea BlkBassStd350_Subarea.weightedYear <- BlkBassStd350_FinalModels_2002 %>% 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 = 'Annual Wt. Avg.') # Summarize Results Table_4.10 <- BlkBassStd350_Subarea.weightedYear %>% # mutate(`_` = 'Black Bass Hg Std. 350mm') %>% select(Subarea, Method, Stdz350) %>% rename(`Black Bass Hg Std. 350mm` = Stdz350) %>% # pivot_wider(names_from=c(Subarea), values_from=c(Stdz350)) %>% signif_df(., 3) setwd(wd) writexl::write_xlsx(Table_4.10, 'Reeval_Impl_Goals_Linkage_Analysis/eval1_TLG Eval & BB Impl Goal/Output/Table 4.10.xlsx') Table_4.10 %>% kbl(caption='Reproduce TMDL Staff Report Table 4.10') %>% kable_paper('hover', full_width=F) ``` ### Regression 2: Black Bass Implementation Goal #### 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"} BBSafeConcs_SubareaWtYrModels <- TLG_BlkBass_Subarea.weightedYear %>% group_by(TLG, Predator) %>% group_modify(~ predictionModels(.x, .knownYcol=Stdz350, .knownXcol=CompositeWeightedAvg, knownX=.x$TLGtarget[1], preferExtrapModel=NULL)) %>% rename(ProtectiveBBConc = predictedY) %>% ungroup ``` #### View Black Bass Std. 350mm Predicted Safe Concentration Regressions <b>`r paste(nrow(BBSafeConcs_SubareaWtYrModels), 'Total Regressions Generated ~ distinct(TLG) would select', nrow(BBSafeConcs_SubareaWtYrModels %>% distinct(TLG)), 'regressions (1 for each TLG)' )`</b> ```{r} # Check Suspect Regressions (identified by filter() conditions) BBSafeConcs_SubareaWtYrModels%>% select(TLG, Predator, ModelName, ProtectiveBBConc, 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(BBSafeConcs_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({ BBSafeConcs_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=ProtectiveBBConc, .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 <- BBSafeConcs_SubareaWtYrModels %>% filter(TLG == 'TL3 Fish (50 - <150 mm)' & ModelName == 'exp') # Replace select lowest SER models with hand selected models BBSafeConcs_SubareaWtYr <- BBSafeConcs_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 BBSafeConcs_SubareaWtYr %>% select(-data, -starts_with("known"), -Model) %>% signif_df(., 3) %>% kbl(caption='Final Black Bass Safe Concentration Selected Models') %>% kable_paper('hover', full_width=F) # Avg SER of all Models avgSER_SubareaWtYr <- BBSafeConcs_SubareaWtYr %>% distinct(TLG, SER) %>% summarise(mean(SER)) %>% pull() avgSER_SubareaWtYr # Select Lowest ProtectiveBBConc as Imp. Goal BBSafeConc_WtYr_final <- BBSafeConcs_SubareaWtYr %>% group_by(knownYcol) %>% filter(ProtectiveBBConc == min(ProtectiveBBConc)) ``` ### Save Black Bass Safe COncentration Model Dataframe to create Figs of Regression Models ```{r} save(BBSafeConcs_SubareaWtYr, file=paste0(wd, '/Reeval_Impl_Goals_Linkage_Analysis/eval1_TLG Eval & BB Impl Goal/Output/Regression Models_Predicted BB Protectiv Conc.RData')) ``` ### Final Black Bass Implementation Goal for Linkage Analysis Board staff selected the most protective Black Bass Standard 350mm Hg concentration of <b>`r signif(BBSafeConc_WtYr_final %>% pull(ProtectiveBBConc), 3)`</b> as the Black Bass Implementation Goal to be used in the linkage analysis. #### 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 <- BBSafeConcs_SubareaWtYr %>% select(TLG, Predator, knownX, ProtectiveBBConc) %>% left_join(., TL4targets, by=c('TLG', 'Predator')) %>% rename(Target = knownX) %>% select(TLG, Predator, Target, PredictTL4_150_500, ProtectiveBBConc) %>% 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) setwd(wd) writexl::write_xlsx(Table_4.9, 'Reeval_Impl_Goals_Linkage_Analysis/eval1_TLG Eval & BB Impl Goal/Output/Table 4.9.xlsx') ``` #### Final Black Bass Standard 350mm Hg Implementation Goal ```{r} title <- '' xTitle <- paste("Wt. Avg. Hg Conc. of", BBSafeConc_WtYr_final$TLG[1],'(mg/Kg)') yTitle <- 'Wt. Avg. Hg Conc. of Standardized 350 mm Black Bass (mg/Kg)' # extract data and model dataframe <- BBSafeConc_WtYr_final$data[[1]] dfKnownX <- BBSafeConc_WtYr_final$knownX # Regression Model model <- BBSafeConc_WtYr_final$Model[[1]] modelName <- BBSafeConc_WtYr_final$ModelName model_predict <- tibble(CompositeWeightedAvg = seq(floor( min(dataframe %>% pull(CompositeWeightedAvg), dfKnownX)*95)/100, ceiling(max(dataframe %>% pull(CompositeWeightedAvg), dfKnownX)*105)/100, length.out=500)) model_predictOrig <- model_predict %>% mutate(Stdz350 = case_when(modelName == 'pwr' | modelName == 'exp' ~ exp(predict(model, model_predict %>% select(CompositeWeightedAvg))), T ~ predict(model, model_predict %>% select(CompositeWeightedAvg)))) #ggplot shapes & colors Subareas <- sort(unique(fishData$Subarea)) mySubareashapes <- c(15, 16, 17, 18, 8, 22, 21, 23, 24, 4) Subareashapes <- setNames(mySubareashapes[1:length(Subareas)], Subareas) # color safe palette - https://personal.sron.nl/~pault/#sec:qualitative mySubareacolors <- c('#0077BB','#33BBEE','#009988','#EE7733','#CC3311','#EE3377', '#919191') # Changed '#BBBBBB' light gray to '#919191' darker gray Subareacolors <- setNames(mySubareacolors[1:length(Subareas)], Subareas) used_colors_orig <- Subareacolors[unique(dataframe %>% pull(Subarea))] used_shapes_orig <- Subareashapes[unique(dataframe %>% pull(Subarea))] # Graph Limits xMax <- ceiling(max(model_predictOrig[[1]])*10)/10 yMax <- ceiling(max(model_predictOrig[[2]], max(dataframe %>% select(Stdz350)))*10)/10 #Plot Data Points & Regression Line regressionGraph <- ggplot() + ggtitle(title) + xlab(xTitle) + ylab(yTitle) + #Data Points geom_point(data=dataframe, aes(x=CompositeWeightedAvg, y=Stdz350, color=Subarea, shape=Subarea), alpha=1, size=4, stroke=1.5, show.legend=T) + #Regression Line geom_point(data=model_predictOrig, aes(x=CompositeWeightedAvg, y=Stdz350), color='grey26', size=.3, show.legend=F) + # Regression Equation annotate("text", x=.15, y=.9, label=paste0("italic(y)==", round(exp(model$coefficients['(Intercept)']),3), "~'*'~", round(exp(model$coefficients['CompositeWeightedAvg']),3), " ^ italic(X)"), parse=T) + #Regression Model SER (Standard Error of Regression) annotate("text", x=.15, y=.84, # BBSafeConc_WtYr_final %>% pull(ProtectiveBBConc)*1.9, label=paste0("italic(~SER) == ", sigfig(BBSafeConc_WtYr_final$SER,3)), parse=T) + #Scale & Theme scale_x_continuous(expand = c(0, 0), limits = c(0, xMax)) + scale_y_continuous(expand = c(0, 0), limits = c(0, yMax)) + scale_color_manual(values = used_colors_orig) + scale_shape_manual(values = used_shapes_orig) + theme_light() + theme(text = element_text(size=14, family="sans"), axis.title=element_text(size=12, face="bold"), axis.text.x = element_text(size=10.5, angle = 0, hjust = 0.5), panel.grid.major.y=element_line(size=1), panel.grid.minor=element_blank(), legend.text = element_text(size = 14), legend.title = element_text(size = 14, face="bold")) regressionGraph #Add Prediction Line predictionGraph <- regressionGraph + #Vertical Line annotate(geom="segment", size=0.8, arrow=arrow(type="open", length=unit(0.02, "npc")), x=BBSafeConc_WtYr_final$knownX, xend=BBSafeConc_WtYr_final$knownX, y=0, yend=BBSafeConc_WtYr_final %>% pull(ProtectiveBBConc), color='grey42') + #Vertical Line Text - Known X annotate("text", x=0.085, hjust=0, y=0.03, label=paste0("italic(x)== ", BBSafeConc_WtYr_final$knownX, "[~",BBSafeConc_WtYr_final$Predator," ~TLG~Target]"), parse=T, color='grey42') + #Horizontal Line annotate(geom="segment", size=0.8, arrow=arrow(type="open", length=unit(0.02, "npc")), x=BBSafeConc_WtYr_final$knownX, xend=0, y=BBSafeConc_WtYr_final %>% pull(ProtectiveBBConc), yend=BBSafeConc_WtYr_final %>% pull(ProtectiveBBConc), color='grey42') + #Horizontal Line Text - Predicted Y annotate("text", x=.04, y=.29, label=paste0('italic(y)== ', sigfig(BBSafeConc_WtYr_final %>% pull(ProtectiveBBConc), 3)), parse=T, color='grey42') predictionGraph ggsave(filename = paste0(wd, '/Reeval_Impl_Goals_Linkage_Analysis/eval1_TLG Eval & BB Impl Goal/Output/Fig 4.4 Black Bass Implementation Goal.png'), plot = predictionGraph, width = 9, height = 5.75, units = "in", dpi = 125) ``` #### Save BB Data Used in final Eval and Implementation Goal for Linkage Analysis ```{r} setwd(wd) # Save Final Black Bass Hg Implementation Goal BBImpGoal <- signif(BBSafeConc_WtYr_final %>% pull(ProtectiveBBConc), 3) save(BBImpGoal, file="Reeval_Impl_Goals_Linkage_Analysis/eval1_TLG Eval & BB Impl Goal/Output/BBImpGoal.RData") writexl::write_xlsx(tibble('BBImpGoal'=BBImpGoal), 'Reeval_Impl_Goals_Linkage_Analysis/eval1_TLG Eval & BB Impl Goal/Output/BBImpGoal.xlsx') ``` ## Save Data used for Final Evaluations ```{r} TLG_save <- TLG %>% 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) # concatenate and save data used for final black bass models in BB Eval bb_eval_save <- NULL for(i in 1:nrow(BlkBassStd350_FinalModels_2002)){ temp <- BlkBassStd350_FinalModels_2002$data[[i]] %>% select(SourceID, ProjectCode, StationName, SampleDate, Analyte, Result, MDL, RL, ResultQualCode, Censored, CommonName, TLAvgLength, TissueName, NumberFishPerComp) %>% rename(ReferenceCode=ProjectCode, TLAvgLength_mm = TLAvgLength) %>% mutate(Subarea = BlkBassStd350_FinalModels_2002$Subarea[i], .before = SourceID) %>% mutate(Unit = BlkBassStd350_FinalModels_2002$Unit[i], .after = Analyte) %>% mutate(SourceID = case_when(SourceID == "CEDENFish" ~ "CEDEN", T ~ SourceID), ReferenceCode = case_when(SourceID == "R5MF" ~ ReferenceCode, T ~ NA_character_)) %>% arrange(Subarea, SampleDate, StationName) bb_eval_save <- bind_rows(bb_eval_save, temp) } writexl::write_xlsx(list("TLG Final Eval"=TLG_save, "Black Bass Final Eval"=bb_eval_save), paste0(wd, '/Reeval_Impl_Goals_Linkage_Analysis/eval1_TLG Eval & BB Impl Goal/Output/Appdx A_Data Compilation_Data Used for Final 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:49:38 PST"