In noaa-garfo/discaRd: Cochran bycatch estimation

knitr::opts_chunk$set(echo=FALSE, warning = FALSE, 
                                            message = FALSE, cache = FALSE,
                                            progress = TRUE, verbose = FALSE, comment = F
                                            , error = FALSE, dev = 'png', dpi = 200)

library(odbc)
library(dplyr, warn.conflicts = FALSE)
# library(dbplyr)
library(ggplot2)
# library(config)
library(stringr)
library(discaRd)
options(scipen = 999)

source('cams_discard_functions.R')

# local run
dw_apsd <- config::get(value = "apsd", file = "K:/R_DEV/config.yml")

# if on server..
#dw_apsd <- config::get(value = "apsd", file = "~/config.yml")

bcon <- dbConnect(odbc::odbc(), 
                                    DSN = dw_apsd$dsn, 
                                    UID = dw_apsd$uid, 
                                    PWD = dw_apsd$pwd)

'%!in%' <- function(x,y)!('%in%'(x,y))

species_nespp3 = '352'  # spiny dogfish

# define species stock if needed
# species_stock = 'GOM'  # GOM cod

# if using a unit stock, make this NULL!!
species_stock = NA  # all unit stocks

FY <- 2018

# Stat areas table

CAMS_STATAREA_STOCK = tbl(bcon, sql('select * from apsd.CAMS_STATAREA_STOCK')) 
# %>% collect()

# GEAR TABLE

  CAMS_GEAR_STRATA = tbl(bcon, sql('  select *   from APSD.CAMS_GEARCODE_STRATA')) %>% 
    collect()

# Mortality table
CAMS_DISCARD_MORTALITY_STOCK = tbl(bcon, sql("select *
from apsd.CAMS_DISCARD_MORTALITY_STOCK")) 
# %>% 
#   collect()

# unique stat areas for stock ID if needed
STOCK_AREAS = CAMS_STATAREA_STOCK %>% 
  filter(NESPP3 == species_nespp3) %>%  # removed  & AREA_NAME == species_stock
  distinct(STAT_AREA) %>% 
  collect()

CAMS_DISCARD_MORTALITY_STOCK = CAMS_DISCARD_MORTALITY_STOCK %>% 
  collect() %>% 
  mutate(SPECIES_STOCK = AREA_NAME) %>% 
  select(-AREA_NAME) %>% 
  filter(NESPP3 == species_nespp3) ## don't reallt want this here...

# get catch and matched obs data together

c_o_dat2 <- tbl(bcon, sql(paste0("
with obs_cams as (
   select year
    , month
    --, case when region = 'N' then 'NE'
  --       when region = 'S' then 'MA'
  --       end as region
    , case when month in (5,6,7,8,9,10) then 1
           when month in (11,12,1,2,3,4) then 2
           end as halfofyear
    , carea
    , vtrserno
    , link1
    , docid
    , dmis_trip_id
    , nespp3
    , GEARCODE
    , NEGEAR
    , GEARTYPE
    , MESHGROUP
    , SECTOR_ID
    , tripcategory
    , accessarea
    , activity_code_1
    , permit_EFP_1
  , permit_EFP_2
  , permit_EFP_3
  , permit_EFP_4
  , redfish_exemption
    , closed_area_exemption
    , sne_smallmesh_exemption
    , xlrg_gillnet_exemption
    , NVL(sum(discard),0) as discard
    , NVL(round(max(subtrip_kall)),0) as subtrip_kall
    , NVL(round(max(obs_kall)),0) as obs_kall
    , NVL(sum(discard)/round(max(obs_kall)), 0) as dk
    from apsd.bg_cams_obs_catch
--  where nespp3 is not null
    group by year, carea, vtrserno, link1, nespp3, docid, GEARCODE, NEGEAR, GEARTYPE
    , MESHGROUP, dmis_trip_id, month
    --, region
    , halfofyear
    , sector_id
    , tripcategory
    , accessarea
    , activity_code_1
    , permit_EFP_1
  , permit_EFP_2
  , permit_EFP_3
  , permit_EFP_4
  , redfish_exemption
    , closed_area_exemption
    , sne_smallmesh_exemption
    , xlrg_gillnet_exemption
    order by vtrserno asc
    ) 

, regions as (

select STAT_AREA 
    , AREA_NAME as REGION
    from APSD.CAMS_STATAREA_STOCK
    where NESPP3 = ", species_nespp3, "
)

, strata as ( 

select CAMS_GEAR_GROUP
     , VTR_GEAR_CODE
    from APSD.CAMS_GEARCODE_STRATA
    where NESPP3 = ", species_nespp3, "
)
--
--select count(distinct(vtrserno)) n_vtr
--, count(distinct(DMIS_TRIP_ID)) n_dmisid
--, count(distinct(LINK1)) n_link1
--, count(distinct(STRATA)) n_strata
----select distinct(GEARCODE)
--from (

     select case when MONTH in (1,2,3,4) then YEAR-1 else YEAR end as GF_YEAR
  , case when MONTH in (1,2,3) then YEAR-1 else YEAR end as SCAL_YEAR
  , o.*
  , r.region
  , c.match_nespp3
  , coalesce(c.match_nespp3, o.nespp3) as nespp3_final
  , NVL(s.CAMS_GEAR_GROUP, '0')||'_'||o.MESHGROUP||'_'||NVL(r.REGION,'na')||'_'||o.HALFOFYEAR as STRATA
  , NVL(s.CAMS_GEAR_GROUP, '0') CAMS_GEAR_GROUP
  from obs_cams o
  left join apsd.s_nespp3_match_conv c on o.nespp3 = c.nespp3
  left join (select * from strata) s 
  ON s.VTR_GEAR_CODE = o.GEARCODE
  left join (select * from regions) r 
  ON r.STAT_AREA = o.CAREA

"
    )
  )
) 
# %>% 
#   collect()

#---------------------------------------------------------------------------------------#
# Here is where various stratification would make sense
# DOCID is hard coded in discaRd.. in this case, CV is calculated using  N = subtrips
#---------------------------------------------------------------------------------------#

ddat_focal <- c_o_dat2 %>% 
  filter(GF_YEAR == FY) %>%   ## time element is here!!
  filter(CAREA %in% local(STOCK_AREAS$STAT_AREA)) %>% 
  # group_by(STRATA) %>% 
  collect() %>% 
  mutate(LIVE_POUNDS = SUBTRIP_KALL
         ,SEADAYS = 0
             , NESPP3 = NESPP3_FINAL) %>% 
  left_join(., y = CAMS_DISCARD_MORTALITY_STOCK %>% 
              dplyr::select(-NESPP3) %>% 
              mutate(REGION = SPECIES_STOCK)
            , by = c('REGION', 'CAMS_GEAR_GROUP')
            ) %>% 
  relocate('COMMON_NAME','SPECIES_ITIS','NESPP3','REGION','CAMS_GEAR_GROUP','Discard_Mortality_Ratio')


ddat_prev <- c_o_dat2 %>% 
  filter(GF_YEAR == FY-1) %>%    ## time element is here!!
  filter(CAREA %in% local(STOCK_AREAS$STAT_AREA)) %>% 
  # group_by(STRATA) %>% 
  collect() %>% 
  mutate(LIVE_POUNDS = SUBTRIP_KALL
         ,SEADAYS = 0
             , NESPP3 = NESPP3_FINAL) %>% 
  left_join(., y = CAMS_DISCARD_MORTALITY_STOCK %>% 
              dplyr::select(-NESPP3) %>% 
              mutate(REGION = SPECIES_STOCK)
            , by = c('REGION', 'CAMS_GEAR_GROUP')
            ) %>% 
  relocate('COMMON_NAME','SPECIES_ITIS','NESPP3','REGION','CAMS_GEAR_GROUP','Discard_Mortality_Ratio')

save(ddat_focal,file=paste0("ddat_focal_nespp",species_nespp3,"_FY",FY,".Rdata"))
save(ddat_prev,file=paste0("ddat_prev_nespp",species_nespp3,"_FY",FY,".Rdata"))

For species that have the same stock definition (i.e. unit stocks), these data pieces may be pulled at the same time. The next set of chunks can be run in sequence on a series of nespp3 (species_nespp3) codes.

load(file=paste0("ddat_focal_nespp",species_nespp3,"_FY",FY,".Rdata"))
load(file=paste0("ddat_prev_nespp",species_nespp3,"_FY",FY,".Rdata"))

# FIX THIS!!! this could be the duping KALL issue.. 

# need to slice the first record for each observed trip.. these trips are multi rowed while unobs trips are single row.. 
ddat_focal_spp = ddat_focal %>% 
  filter(!is.na(LINK1)) %>% 
  group_by(LINK1, VTRSERNO) %>% 
  slice(1) %>% 
  ungroup()

# and join to the unobserved trips

ddat_focal_spp = ddat_focal_spp %>% 
  union_all(ddat_focal %>% 
              filter(is.na(LINK1)) %>% 
               group_by(VTRSERNO) %>% 
               slice(1) %>% 
               ungroup()
            )


# if using the combined catch/obs table, which seems necessary for groundfish.. need to roll your own table to use with run_discard function

bdat_spp = ddat_focal %>% 
  filter(!is.na(LINK1)) %>% 
  mutate(DISCARD_PRORATE = DISCARD
         , OBS_AREA = CAREA
         , OBS_HAUL_KALL_TRIP = OBS_KALL
         , PRORATE = 1)

ddat_discards = run_discard(bdat = bdat_spp
                                             , ddat = ddat_focal_spp
                                             , c_o_tab = ddat_focal
                                             # , year = 2019
                                             # , species_nespp3 = '081' # haddock...
                                             , species_nespp3 = species_nespp3  #'081' #cod... 
                                             , stratvars = c('CAMS_GEAR_GROUP', 'MESHGROUP', 'REGION', 'HALFOFYEAR')
                                              # , stratvars = c('GEARTYPE', 'MESHGROUP', 'REGION', 'HALFOFYEAR')  #OLD GEAR TYPE
                                             , aidx = c(1,2)
                                             )

ddat_discards$res$Discard_Mortality_Ratio[is.na(ddat_discards$res$Discard_Mortality_Ratio)] <- 1

# sum of all discard, by subtrip, multiplied by mortality ratio
sum(ddat_discards$res$DISCARD*ddat_discards$res$Discard_Mortality_Ratio, na.rm = T) # 

# discard rates by strata
dest_strata = ddat_discards$allest$C %>% summarise(STRATA = STRATA
                       , N = N
                       , n = n
                       , orate = round(n/N, 2)
                       , drate = RE_mean
                       , KALL = K, disc_est = round(D)
                       , CV = round(RE_rse, 2)
                                        )

# rates for trawls
dest_strata %>% slice(grep('50_*', dest_strata$STRATA))

ddat_discards$res %>%  mutate(newd = DISCARD*Discard_Mortality_Ratio) %>% group_by(STRATA) %>% dplyr::summarise(D=sum(newd,na.rm=T),d_rate=mean(DISC_RATE,na.rm=T)) %>% arrange(desc(D)) %>% print.data.frame()

# Potentially limited to the focal year (for year end estimation)

# create a small bdat
bdat_mod <- bdat_spp %>% filter(
  OBS_KALL > 0,
  NESPP3_FINAL==species_nespp3
  ) %>%
  select(DMIS_TRIP_ID,SUBTRIP_KALL,OBS_KALL,DISCARD,
         STRATA,CAMS_GEAR_GROUP,REGION,HALFOFYEAR,MESHGROUP)

# full stratification (similar to ratio estimator)
mod.full <- glm(DISCARD~offset(log(OBS_KALL)) + STRATA -1,
           family=poisson(),data=bdat_mod)
# additive effects (to accommodate unobserved full-factorial strata)
mod <- glm(DISCARD~offset(log(OBS_KALL)) + CAMS_GEAR_GROUP + MESHGROUP + HALFOFYEAR + REGION,
           family=poisson(),data=bdat_mod)

# create ddat for model predictions
ddat_pred <- ddat_discards$res %>% 
  filter(
    #STRATA %in% bdat_spp_mod$STRATA,  #only for full factorial
    SUBTRIP_KALL > 0
    ) %>%
  select(DMIS_TRIP_ID,SUBTRIP_KALL,STRATA,Discard_Mortality_Ratio,EST_DISCARD,
         CAMS_GEAR_GROUP,REGION,HALFOFYEAR,MESHGROUP) %>%
  mutate(
    OBS_KALL = SUBTRIP_KALL,
    # adjust for unobserved gear groups
    CAMS_GEAR_GROUP =
           case_when(
             CAMS_GEAR_GROUP %in% unique(bdat_mod$CAMS_GEAR_GROUP) ~ CAMS_GEAR_GROUP,
             TRUE ~ '0'
           ),
    NULL=NULL) 

# model predictions
ddat_pred <- data.frame(ddat_pred,D_MODEL=predict(mod,ddat_pred,type = "response"))
# sum predicted discards for all trips
sum(ddat_pred$D_MODEL)

# sum of predicted discards for:
#    1) trips with unobserved strata (model)
#    2) trips with observed strata (ratio estimator)

ddat_pred %>%
  mutate(DEAD_DISCARDS_FINAL =
           case_when(
             STRATA %in% dest_strata$STRATA[is.na(dest_strata$drate)] ~ D_MODEL*Discard_Mortality_Ratio,
             STRATA %!in% dest_strata$STRATA[is.na(dest_strata$drate)] ~ EST_DISCARD*Discard_Mortality_Ratio
           )) %>%
  group_by(REGION) %>%
  dplyr::summarise(sum(DEAD_DISCARDS_FINAL))

The next chunk could be combined with the previous one if running several species in sequence.

# set up trips table for current year
# done above 


# set up trips table for previous year
ddat_prev_spp = ddat_prev %>% 
  filter(!is.na(LINK1)) %>% 
  group_by(LINK1) %>% 
  slice(1)

ddat_prev_spp = ddat_prev_spp %>% 
  union_all(ddat_prev %>% filter(is.na(LINK1)))


# previous year observer data needed.. 
bdat_prev_spp = ddat_prev %>% 
  filter(!is.na(LINK1)) %>% 
  mutate(DISCARD_PRORATE = DISCARD
         , OBS_AREA = CAREA
         , OBS_HAUL_KALL_TRIP = OBS_KALL
         , PRORATE = 1)

# Run the discaRd functions on previous year
d_prev = run_discard(bdat = bdat_prev_spp
                                             , ddat = ddat_prev_spp
                                             , c_o_tab = ddat_prev
                                             # , year = 2018
                                             , species_nespp3 = species_nespp3
                                             , stratvars = c('CAMS_GEAR_GROUP', 'MESHGROUP', 'REGION', 'HALFOFYEAR')
                                             , aidx = c(1,2)
                                             )


# Run the discaRd functions on current year
d_focal = run_discard(bdat = bdat_spp
                                             , ddat = ddat_focal_spp
                                             , c_o_tab = ddat_focal
                                             # , year = 2019
                                             # , species_nespp3 = '081' # haddock...
                                             , species_nespp3 = species_nespp3  #'081' #cod... 
                                             , stratvars = c('CAMS_GEAR_GROUP', 'MESHGROUP', 'REGION', 'HALFOFYEAR')
                                             , aidx = c(1,2)
                                             )

# summarize each result for convenience
dest_strata_p = d_prev$allest$C %>% summarise(STRATA = STRATA
                       , N = N
                       , n = n
                       , orate = round(n/N, 2)
                       , drate = RE_mean
                       , KALL = K, disc_est = round(D)
                       , CV = round(RE_rse, 2)
                                        )

dest_strata_f = d_focal$allest$C %>% summarise(STRATA = STRATA
                       , N = N
                       , n = n
                       , orate = round(n/N, 2)
                       , drate = RE_mean
                       , KALL = K, disc_est = round(D)
                       , CV = round(RE_rse, 2)
                                        )

# substitute transition rates where needed

trans_rate_df = dest_strata_f %>% 
  left_join(., dest_strata_p, by = 'STRATA') %>% 
  mutate(STRATA = STRATA
         , n_obs_trips = n.x
         , in_season_rate = drate.x
         , previous_season_rate = drate.y
         , trans_rate = get.trans.rate(l_observed_trips = n_obs_trips
                                         , l_assumed_rate = previous_season_rate
                                         , l_inseason_rate = in_season_rate
                                         )
         ) %>% 
  dplyr::select(STRATA
         , n_obs_trips 
         , in_season_rate 
         , previous_season_rate 
         , trans_rate)


trans_rate_df = trans_rate_df %>% 
  mutate(final_rate = case_when((in_season_rate != trans_rate & !is.na(trans_rate)) ~ trans_rate)) 

 trans_rate_df$final_rate = coalesce(trans_rate_df$final_rate, trans_rate_df$in_season_rate)

# compare subbed rate total to in season total. 
 trans_rate_df %>% 
   left_join(dest_strata_f, by = 'STRATA') %>% 
   summarise(in_season_discard = sum(in_season_rate*KALL, na.rm = T)
             , final_rate_discard = sum(final_rate*KALL, na.rm = T)
             ) %>% 
   knitr::kable()

# park the final rate into the trips table
 trans_rate_df %>% 
   right_join(., y = d_focal$res, by = 'STRATA') %>% 
   summarise(inseason_rate_d = sum(in_season_rate*SUBTRIP_KALL*Discard_Mortality_Ratio, na.rm = T)
             , trans_rate_d = sum(trans_rate*SUBTRIP_KALL*Discard_Mortality_Ratio, na.rm = T)
             , final_rate_d = sum(final_rate*SUBTRIP_KALL*Discard_Mortality_Ratio, na.rm = T))

assumed_discard %>% knitr::kable(caption = 'Assumed discard. This could be a generalized rate from current year. This could also be built from previous years discard info.')

dest_strata %>% knitr::kable(caption = 'Stratified Estimate From discaRd')

# out_tab %>% head() %>% 
    # knitr::kable(format= 'markdown', caption = "Exampe of tabular output")