Introduction

Here we briefly test changes to the function calc_age2length that should allow the user to specify a different max length bin. For all setup, etc, please see previous files Full methods are explained here and here.

This page has visualizations for the NOBA model example, CERES Global Sustainability. For full explanation of methods, see the file linked at the beginning of each section.

Simulate a survey part 4: sample for length composition (testing revised function)

Full methods are explained here.

We will apply examples here to only one species, Greenland halibut in NOBA, which grows to a large size.

To create a survey, the user specifies the timing of the survey, which species are captured, the spatial coverage of the survey, the species-specific survey efficiency (“q”), and the selectivity at age for each species.

# make defaults that return a standard survey, implement in standard_survey
# users need to map their species groups into these general ones
#   large pelagics/reef associated/burrowers/otherwise non-trawlable
#   pelagics
#   demersals
#   selected flatfish
if(initNOBA) funct.groups <- rename(funct.groups, GroupType = InvertType)

survspp <- funct.groups$Name[funct.groups$IsTurnedOn==1 &
                           funct.groups$GroupType %in% c("FISH", "SHARK")]

if(initCCA) { #Sarah's CCA Grouping
  nontrawl <- c("Shark_C","Yelloweye_rockfish","Benthopel_Fish","Pisciv_S_Fish",
                "Pisciv_T_Fish","Shark_D","Shark_P")
  pelagics <- c("Pisciv_V_Fish","Demersal_S_Fish","Pacific_Ocean_Perch","Mesopel_M_Fish",
                "Planktiv_L_Fish","Jack_mackerel","Planktiv_S_Fish","Pacific_sardine",
                "Anchovy","Herring","Pisciv_B_Fish")
  demersals <- c("Demersal_P_Fish","Planktiv_O_Fish","Demersal_D_Fish",
                 "Demersal_DC_Fish","Demersal_O_Fish","Darkblotched_rockfish",
                 "Demersal_F_Fish","Demersal_E_Fish","Bocaccio_rockfish",
                 "Demersal_B_Fish","Shark_R","Mesopel_N_Fish","Shark_B","Spiny_dogfish",
                 "SkateRay")
  selflats <- c("Pisciv_D_Fish", "Arrowtooth_flounder","Petrale_sole")
}

if(initNEUS) { # Sarah's NEUS Grouping
  nontrawl <- c("Pisciv_T_Fish", "Shark_D", "Shark_P", "Reptile", "Mesopel_M_Fish")
  pelagics <- c("Planktiv_L_Fish", "Planktiv_S_Fish", "Benthopel_Fish", "Pisciv_S_Fish")
  demersals <- c("Pisciv_D_Fish", "Demersal_D_Fish","Demersal_E_Fish", 
                 "Demersal_S_Fish","Demersal_B_Fish","Demersal_DC_Fish",
                 "Demersal_O_Fish","Demersal_F_Fish",
                 "Shark_B", "SkateRay")
  selflats <- c("Pisciv_B_Fish")
}

if(initNOBA) { # Sarah's NOBA Grouping
  nontrawl <- c("Sharks_other", "Pelagic_large","Mesop_fish")
  pelagics <- c("Pelagic_small","Redfish_other","Mackerel","Haddock",
                "Saithe","Redfish","Blue_whiting","Norwegian_ssh","Capelin")
  demersals <- c("Demersals_other","Demersal_large","Flatfish_other","Skates_rays",
                 "Green_halibut","North_atl_cod","Polar_cod","Snow_crab")
  selflats <- c("Long_rough_dab")
}

The following settings are for our example standard survey once per year, most areas, with mixed efficiency and selectivity:

# general specifications for bottom trawl survey, with items defined above commented out to avoid wasting time loading already loaded files:
#   once per year at mid year
# generalized timesteps all models
runpar <- load_runprm(d.name, run.prm.file)
noutsteps <- runpar$tstop/runpar$outputstep
stepperyr <- if(runpar$outputstepunit=="days") 365/runpar$toutinc

midptyr <- round(median(seq(1,stepperyr)))

annualmidyear <- seq(midptyr, noutsteps, stepperyr)

#   ~75-80% of boxes (leave off deeper boxes?)
boxpars <- load_box(d.name, box.file)
boxsurv <- c(2:round(0.8*(boxpars$nbox - 1)))

#   define bottom trawl mixed efficiency
ef.nt <- 0.01 # for large pelagics, reef dwellers, others not in trawlable habitat
ef.pl <- 0.1  # for pelagics
ef.dm <- 0.7  # for demersals
ef.fl <- 1.1  # for selected flatfish

# bottom trawl survey efficiency specification by species group
effnontrawl <- data.frame(species=nontrawl, efficiency=rep(ef.nt,length(nontrawl)))
effpelagics <- data.frame(species=pelagics, efficiency=rep(ef.pl,length(pelagics)))
effdemersals <- data.frame(species=demersals, efficiency=rep(ef.dm,length(demersals)))
effselflats <- data.frame(species=selflats, efficiency=rep(ef.fl,length(selflats)))

efficmix <- bind_rows(effnontrawl, effpelagics, effdemersals, effselflats)

#   mixed selectivity (using 10 agecl for all species)
#     flat=1 for large pelagics, reef dwellers, others not in trawlable habitat
#     sigmoid 0 to 1 with 0.5 inflection at agecl 3 for pelagics, reaching 1 at agecl 5, flat top
#     sigmoid 0 to 1 with 0.5 inflection at agecl 5 for most demersals and flatfish, reaching 1 at agecl 7, flat top
#     dome shaped 0 to 1 at agecl 6&7 for selected demersals, falling off to 0.7 by agecl 10

sigmoid <- function(a,b,x) {
  1 / (1 + exp(-a-b*x))
}

# survey selectivity specification by species group
selnontrawl <- data.frame(species=rep(nontrawl, each=10),
                          agecl=rep(c(1:10),length(nontrawl)),
                          selex=rep(1.0,length(nontrawl)*10))
selpelagics <- data.frame(species=rep(pelagics, each=10),
                          agecl=rep(c(1:10),length(pelagics)),
                          selex=sigmoid(5,1,seq(-10,10,length.out=10)))
seldemersals <- data.frame(species=rep(demersals, each=10),
                          agecl=rep(c(1:10),length(demersals)),
                          selex=sigmoid(1,1,seq(-10,10,length.out=10)))
selselflats <- data.frame(species=rep(selflats, each=10),
                          agecl=rep(c(1:10),length(selflats)),
                          selex=sigmoid(1,1,seq(-10,10,length.out=10)))

selexmix <- bind_rows(selnontrawl, selpelagics, seldemersals, selselflats)

# use this constant 0 cv for testing
surv_cv_0 <- data.frame(species=survspp, cv=rep(0.0,length(survspp)))

#   define bottom trawl survey cv by group
cv.nt <- 1.0 # for large pelagics, reef dwellers, others not in trawlable habitat
cv.pl <- 0.5  # for pelagics
cv.dm <- 0.3  # for demersals
cv.fl <- 0.3  # for selected flatfish

# specify cv by species groups
surv_cv_nontrawl <- data.frame(species=nontrawl, cv=rep(cv.nt,length(nontrawl)))
surv_cv_pelagics <- data.frame(species=pelagics, cv=rep(cv.pl,length(pelagics)))
surv_cv_demersals <- data.frame(species=demersals, cv=rep(cv.dm,length(demersals)))
surv_cv_selflats <- data.frame(species=selflats, cv=rep(cv.fl,length(selflats)))

surv_cv_mix <- bind_rows(surv_cv_nontrawl, surv_cv_pelagics, surv_cv_demersals, surv_cv_selflats)

And the numbers to be sampled for lengths each year:

#   define n fish for biological sampling by group
#   this could easily be a vector or time series, constant here
ns.nt <- 25 # for large pelagics, reef dwellers, others not in trawlable habitat
ns.pl <- 1000  # for pelagics
ns.dm <- 1000  # for demersals
ns.fl <- 1000  # for selected flatfish

effNnontrawl <- data.frame(species=nontrawl, effN=rep(ns.nt,length(nontrawl)))
effNpelagics <- data.frame(species=pelagics, effN=rep(ns.pl,length(pelagics)))
effNdemersals <- data.frame(species=demersals, effN=rep(ns.dm,length(demersals)))
effNselflats <- data.frame(species=selflats, effN=rep(ns.fl,length(selflats)))

effNmix <- bind_rows(effNnontrawl, effNpelagics, effNdemersals, effNselflats)

Here we use create_survey on the numbers output of run_truth to create the survey census of age composition (for just one species in this case). The sample_fish applies the median for aggregation and does not apply multinomial sampling if sample=FALSE in the function call.

ss.name <- funct.group.names[funct.group.names == "Green_halibut"]

#change back to flat selectivity to see full comp
selex1 <- data.frame(species=rep(funct.group.names, each=10),
                     agecl=rep(c(1:10),length(funct.group.names)),
                     selex=rep(1.0,length(funct.group.names)*10))

# get survey nums with full (no) selectivity
ss_survey_testNstd_nosel <- create_survey(dat = truth$nums,
                                 time = annualmidyear,
                                 species = ss.name,
                                 boxes = boxsurv,
                                 effic = efficmix,
                                 selex = selex1)

# now sample fish nums from this
ss_numsstd_nosel <- sample_fish(ss_survey_testNstd_nosel, effNmix)

# aggregate true resn per survey design
survey_aggresnstd <- aggregateDensityData(dat = truth$resn,
                                 time = annualmidyear,
                                 species = ss.name,
                                 boxes = boxsurv)

# aggregate true structn per survey design
survey_aggstructnstd <- aggregateDensityData(dat = truth$structn,
                                 time = annualmidyear,
                                 species = ss.name,
                                 boxes = boxsurv)

#dont sample these, just aggregate them using median (effNmix does nothing)
ss_structnstd <- sample_fish(survey_aggstructnstd, effNmix, sample = FALSE)

ss_resnstd <-  sample_fish(survey_aggresnstd, effNmix, sample = FALSE)

Length sample with default max length bin (150 cm):

ss_length_stdsurv_nosel <- calc_age2length(structn = ss_structnstd,
                                 resn = ss_resnstd,
                                 nums = ss_numsstd_nosel,
                                 biolprm = truth$biolprm, fgs = truth$fgs,
                                 CVlenage = 0.1, remove.zeroes=TRUE)

Length sample with user specified max length bin (250 cm):

ss_length_stdsurv_nosel_max <- calc_age2length(structn = ss_structnstd,
                                 resn = ss_resnstd,
                                 nums = ss_numsstd_nosel,
                                 biolprm = truth$biolprm, fgs = truth$fgs,
                                 maxbin = 250,
                                 CVlenage = 0.1, remove.zeroes=TRUE)

Plots show that default length bins are not adequate for something like Greenland halibut, which get larger than the default largest bin in calc_age2length, 150 cm. Both census and these sampled length comps are chopped off at 150 cm for this species:

I changed the calc_age2length code to pass the upper bin of 150 as a default but allow the user to change it in the function call. This plot shows results from setting maxbin = 250 for Greenland Halibut in NOBA:

Now we can go back and apply standard survey selectivity as in the NOBA visualization and see how the length comps look:

ss.name <- funct.group.names[funct.group.names == "Green_halibut"]

# get survey nums with standard selectivity
ss_survey_testNstd <- create_survey(dat = truth$nums,
                                 time = annualmidyear,
                                 species = ss.name,
                                 boxes = boxsurv,
                                 effic = efficmix,
                                 selex = selexmix)

# now sample fish nums from this
ss_numsstd <- sample_fish(ss_survey_testNstd, effNmix)

# structn and resn stuff is exactly the same because selectivity is irrelevant

Length sample from standard survey with user specified max length bin (250 cm):

ss_length_stdsurv_max <- calc_age2length(structn = ss_structnstd,
                                 resn = ss_resnstd,
                                 nums = ss_numsstd,
                                 biolprm = truth$biolprm, fgs = truth$fgs,
                                 maxbin = 250,
                                 CVlenage = 0.1, remove.zeroes=TRUE)

Bigger max bin is definitely more appropriate for this species.

I suppose we could make the maxbin a vector by species, but for now may be overkill.