Merge branch 'main' of https://github.com/pfmc-assessments/yellowtail…

…_2025

Rscripts/logbook_analys.R

@@ -0,0 +1,87 @@
+library(dplyr)
+
+observer <- read.csv('data/confidential/wcgop_logbook/ward_oken_observer_data_2024-09-10.csv') |>
+  as_tibble() |>
+  filter(AVG_LAT > 40+1/6, 
+         sector == 'Midwater Rockfish',
+         DATATYPE == 'Analysis Data') 
+
+observer_yt <- observer|>
+  group_by(HAUL_ID) |>
+  summarise(has_yt = 'Yellowtail Rockfish' %in% species,
+            across(AVG_DEPTH:YEAR, first)) |>
+  filter(!has_yt) |>
+  mutate(species = 'Yellowtail Rockfish', MT = 0) |>
+  select(-has_yt) |>
+  bind_rows(
+    filter(observer, species == 'Yellowtail Rockfish')
+  ) |>
+  mutate(DRVID = as.integer(DRVID))
+
+em <- read.csv('data/confidential/wcgop_logbook/ward_oken_em_data_2024-12-03.csv') |>
+  as_tibble() |> 
+  filter(AVG_LAT > 40+1/6, sector == 'Midwater Rockfish EM') 
+
+em_yt <- em |>
+  group_by(HAUL_ID) |>
+  summarise(has_yt = 'Yellowtail Rockfish' %in% species,
+            across(AVG_DEPTH:YEAR, first)) |> 
+  filter(!has_yt) |>
+  mutate(species = 'Yellowtail Rockfish', MT = 0) |>
+  select(-has_yt) |>
+  bind_rows(
+    filter(em, species == 'Yellowtail Rockfish')
+  )
+
+all_hauls_yt <- bind_rows(list(obs = observer_yt, em = em_yt), .id = 'program') |>
+  filter(HAUL_DURATION > 0, HAUL_DURATION < 10) |>
+  mutate(depth_std = (AVG_DEPTH - mean(AVG_DEPTH))/sd(AVG_DEPTH),
+         fyear = factor(YEAR), log_duration = log(HAUL_DURATION),
+         cpue = MT/HAUL_DURATION, month = lubridate::month(as.Date(SET_DATE)),
+         fDRVID = factor(DRVID), fport = factor(R_PORT),
+         program = factor(program)) # note D_PORT is NA in some cases, R_PORT is always defined (fish ticket)
+
+nrow(filter(all_hauls_yt, MT>0)) / nrow(all_hauls_yt)
+# good to monitor % positive hauls. This is pretty high (70%).
+
+yt_mesh <- sdmTMB::make_mesh(all_hauls_yt, xy_cols = c('AVG_LONG', 'AVG_LAT'), n_knots = 50)
+plot(yt_mesh)
+# this does not look like a st model is a good idea. very clear "hotspots" of fishing effort.
+# but we need it for sdmTMB to work (not actually used)
+
+mod_simple <- sdmTMB::sdmTMB(formula = cpue ~ fyear + s(depth_std) + s(month, bs = 'cc') + program + (1|fDRVID) + (1|fport) - 1, family = sdmTMB::delta_lognormal(),
+                             data = all_hauls_yt, spatial = 'off', spatiotemporal = 'off', time = 'fyear', mesh = yt_mesh)
+
+summary(mod_simple)
+
+
+resids <- residuals(mod_simple, model = 2) # the presence/absence model looks consistently good
+qqnorm(resids)
+qqline(resids)
+
+sdmTMB::visreg_delta(mod_simple, xvar = "depth_std", model = 1, gg = TRUE)
+sdmTMB::visreg_delta(mod_simple, xvar = "depth_std", model = 2, gg = TRUE)
+sdmTMB::visreg_delta(mod_simple, xvar = "month", model = 1, gg = TRUE)
+sdmTMB::visreg_delta(mod_simple, xvar = "month", model = 2, gg = TRUE)
+
+preds <- predict(mod_simple, 
+                 newdata = tibble(fyear = factor(2012:2023), depth_std = 0, month = 6, program = 'obs', fDRVID = '546053', fport = 'NEWPORT'), 
+                 return_tmb_object = TRUE, re_form = NA, re_form_iid = NA)
+
+ind <- sdmTMB::get_index(preds, bias_correct = TRUE)
+
+ind |> 
+  ggplot(aes(x = as.numeric(fyear), y = est, ymin = lwr, ymax = upr)) +
+  geom_line() +
+  geom_errorbar()
+
+
+all_hauls_yt |>
+  ggplot() +
+  geom_point(aes(x = HAUL_DURATION, y = MT))
+  geom_density(aes(x = HAUL_DURATION, col = YEAR, group = YEAR))
+
+### general observation: as I have added terms to the model, the signal in the index went from well-defined to noisy and meaningless.
+### the port (8 levels) and vessel (46 levels) effects really did it in. should think more carefully about what terms to include.
+### note 3962 hauls, 2762 positive hauls over 12 years. ~330 hauls/yr, 230 positive hauls/yr. # of hauls varies a lot by year.
+### (increases 2012-2018)

docs/data_summary_doc.qmd

@@ -25,37 +25,42 @@ This documents data streams for use in the Northern Yellowtail Rockfish stock as
 Some summary figures:
 
 ```{r pacfin-catch, cache=TRUE}
-load(here('data/confidential/commercial/pacfin_12-09-2024/PacFIN.YTRK.CompFT.09.Dec.2024.RData'))
+load(here('data/confidential/commercial/pacfin_12-11-2024/PacFIN.YTRK.CompFT.11.Dec.2024.RData'))
 # I have checked and LANDING_YEAR = PACFIN_YEAR for all rows.
 
 yt_n_catch_pacfin <- catch.pacfin |>
   filter(
-    AGENCY_CODE == 'W' | AGENCY_CODE == 'O' | COUNTY_NAME == 'DEL NORTE' | COUNTY_NAME == 'HUMBOLDT',
+    AGENCY_CODE == 'W' | AGENCY_CODE == 'O' | COUNTY_NAME == 'DEL NORTE' | COUNTY_NAME == 'HUMBOLDT', # OR, WA, N. CA
     LANDING_YEAR < 2024
   )
 
 # OR landings time series
 or_comm_catch <- read.csv(here("Data/Confidential/commercial/Oregon Commercial landings_433_2023.csv")) |>
   tibble::as_tibble() |>
   group_by(YEAR) |>
-  summarise(catch_mt = sum(TOTAL))
+  summarise(catch_mt = sum(TOTAL)) |>
+  mutate(AGENCY_CODE = 'O')
 
-# WA URCK/UPOP reconstruction
+# WA URCK/UPOP reconstruction, use prior to 2000
 wa_urck <- read.csv(here('Data/Raw_not_confidential/WA_URCK_UPOP_recon_1981-2016.csv')) |>
-  filter(Species_AfterComps %in% c('YTRK', 'YTR1')) |> 
+  filter(Species_AfterComps %in% c('YTRK', 'YTR1'), 
+         Year < 2001) |> 
   group_by(Year) |>
   summarise(catch_mt = sum(RoundPounds_AfterComps) * 0.453592 / 1000) |>
-  mutate(AGENCY_CODE = 'W',
-         model = 'URCK/UPOP')
+  mutate(AGENCY_CODE = 'W') |>
+  rename(YEAR = Year)
 
 # Put it all together
 comm_draft_landings <- yt_n_catch_pacfin |> 
-  filter(AGENCY_CODE != 'O') |> 
+  filter(AGENCY_CODE != 'O',
+         AGENCY_CODE != 'W' | PACFIN_YEAR >= 2001
+         ) |> 
   group_by(PACFIN_YEAR, AGENCY_CODE)|>
   summarise(catch_mt = sum(ROUND_WEIGHT_MTONS)) |>
   rename(YEAR = PACFIN_YEAR) |>
   bind_rows(
-    mutate(or_comm_catch, AGENCY_CODE = 'O')
+    or_comm_catch,
+    wa_urck
   )|> 
   group_by(YEAR, AGENCY_CODE) |>
   summarise(catch_mt = sum(catch_mt)) |>
@@ -109,23 +114,7 @@ comm_draft_landings |>
   geom_line(aes(x = Year, y = catch_mt, col = AGENCY_CODE, linetype = model), linewidth = 1) +
   viridis::scale_color_viridis(discrete = TRUE) +
   labs(title = 'Comparing commercial landings in 2017 model vs current', x = 'YEAR')
-
-# compare 2017 landings to current landings, WA only
-comm_draft_landings |>
-  rename(Year = YEAR) |>
-  mutate(model = 'PacFIN') |>
-  filter(AGENCY_CODE == 'W') |>
-  bind_rows(
-    filter(comm2017, AGENCY_CODE == 'W')
-    )|> 
-  bind_rows(wa_urck) |>
-  ggplot() +
-  geom_line(aes(x = Year, y = catch_mt, col = model), linewidth = 1) +
-  viridis::scale_color_viridis(discrete = TRUE) +
-  labs(title = 'Possible sources for WA landings (vert. line = 1994)', x = 'YEAR') +
-  geom_vline(xintercept = 1994)
 ```
-URCK/UPOP = sum of RoundPounds_AfterComps for Species_AfterComps = YTR1 or YTRK, converted to mt.
 
 Table of commercial landings by state, 1981 onward.
 
@@ -139,9 +128,11 @@ comm_draft_landings |>
   knitr::kable(align = 'l', digits = 1)
 ```
 
-WDFW has alerted the STAT these do not include all tribal catches in recent years.
+WDFW has alerted the STAT these do not include all tribal catches in recent years. 1981-2000 are from the UPOP/URCK reconstruction, summing catches for YTRK and YTR1. 2001 onward comes from PacFIN
 
-Note this is for catch landed into Del Norte and Humboldt counties only.
+Oregon landings were provided directly to the STAT from ODFW and come from PacFIN and a UPOP/URCK reconstruction. 
+
+California landings are for catch landed into Del Norte and Humboldt counties only.
 
 ### At-Sea landings
 
@@ -159,14 +150,15 @@ ashop_catch |>
 
 ```{r rec-catch-modern}
 # WASHINGTON
-wa_modern <- read.csv(here('Data/Raw_not_confidential/RecFIN_WA_catch_to_2023.csv')) |>
-  filter(RECFIN_WATER_AREA_NAME != 'Canada', RECFIN_YEAR <= 2023) |>
+wa_modern <- read.csv(here('Data/Confidential/Rec/CTE501-WASHINGTON-1990---2023.csv')) |> 
+  filter(RECFIN_WATER_AREA_NAME != 'CANADA', RECFIN_YEAR <= 2023, 
+         RECFIN_CATCH_AREA_NAME != 'SEKIU AND PILLAR POINT') |> 
   tibble::as_tibble() |> 
   group_by(RECFIN_YEAR) |> 
-  summarise(Dead_Catch = sum(SUM_TOTAL_MORTALITY_MT)) |>
+  summarise(Dead_Catch = sum(TOTAL_MORTALITY_MT)) |>
   mutate(State = 'WA (mt)')
 
-wa_historical <- read.csv(here('Data/Raw_not_confidential/WA_historical_rec.csv')) |>
+wa_historical <- read.csv(here('Data/Raw_not_confidential/WA_historical_rec.csv')) |> 
   tibble::as_tibble() |>
   filter(AREA %in% 1:4) |> # these are marine areas, per T. Tsou on 12/10/24 in meeting notes
   group_by(RECFIN_YEAR) |>
@@ -247,7 +239,8 @@ ca_mrfss_catch <- read.csv(here('Data/Confidential/Rec/RecFIN_CA_MRFSS.csv')) |>
                              n_ca_catch_mt * weighted_ratio),
          State = 'CA (mt)') |>
   select(-n_ca_catch_mt, -weighted_ratio) |>
-  rename(RECFIN_YEAR = YEAR)
+  rename(RECFIN_YEAR = YEAR) |>
+  filter(RECFIN_YEAR != 1980)
 
 # Interpolate between 3-year average before and after MRFSS break
 ca_mrfss_interpolate <- tibble(
@@ -298,7 +291,7 @@ For CA rec catch estimates for the assessment:
 Initial length sample sizes after running `PacFIN.Utilities::cleanPacFIN()`:
 
 ```{r comm-bio, cache=TRUE}
-load(here('Data/Confidential/Commercial/pacfin_12-09-2024/PacFIN.YTRK.bds.09.Dec.2024.RData'))
+load(here('Data/Confidential/Commercial/pacfin_12-11-2024/PacFIN.YTRK.bds.11.Dec.2024.RData'))
 bds_clean <- PacFIN.Utilities::cleanPacFIN(bds.pacfin, verbose = FALSE) 
 bds_clean |>
   filter(state == 'WA' | state == 'OR' | PACFIN_GROUP_PORT_CODE == 'CCA' |