I am attempting to estimate annual survival by length class (or age class) from the Baranov Catch Equation using r4ss. The specific quantity of interest is exp(-(M + F * S)) where M is natural mortality, F is commercial fishing mortality and S is selectivity. I have written the following R code to perform this operation myself, but am not certain I have done so correctly. I was hoping to compare my approach with the approach used in r4ss. However, I cannot locate where r4ss does this operation. Even if my approach is correct with the example below I am not certain it will be correct if the EXPLOITATION table is creating using different settings.
Below is my R code which calculates annual survival using the LEN_SELEX table and the EXPLOITATION table. I use two approaches. One uses the Dead values from the LEN_SELEX table. The other uses the equation: commercial mortality * selectivity * (retain + (1 - retain) * discard). Both return the same output.
Does r4ss perform the same or a similar operation? If so, could you tell me the function to use? So far I have looked at 'SS_output', 'SS_plots', 'SSbiologytables', 'SSexecutivesummary' and 'SSplotSummaryF'. The function 'SSexecutivesummary' mostly returns errors.
If r4ss does not perform this operation could you identify another R package that does, perhaps ss3sim?
# LEN_SELEX table for 2016
Lsel_table <- read.table(text = '
Factor Fleet Yr Sex Label X49 X51 X53 X55 X57 X59 X61 X63 X65
Lsel 1 2016 1 2016_1_Lsel 0.526523000 0.611541000 0.69631000 0.7772370 0.8505130 0.912405 0.959565 0.989334 0.999987
Lsel 2 2016 1 2016_2_Lsel 0.147894000 0.180769000 0.21850500 0.2611950 0.3087690 0.360966 0.417316 0.477122 0.539458
Ret 1 2016 1 2016_1_Ret 0.061028200 0.138097000 0.28314300 0.4933330 0.7059060 0.855426 0.935835 0.972933 0.988835
Mort 1 2016 1 2016_1_Mort 0.500000000 0.500000000 0.50000000 0.5000000 0.5000000 0.500000 0.500000 0.500000 0.500000
Keep 1 2016 1 2016_1_Keep 0.032132700 0.084451900 0.19715500 0.3834370 0.6003820 0.780495 0.897994 0.962556 0.988822
Dead 1 2016 1 2016_1_Dead 0.279328000 0.347996000 0.44673200 0.5803370 0.7254470 0.846450 0.928780 0.975945 0.994405
Ret 2 2016 1 2016_2_Ret 0.000860749 0.003427850 0.01350010 0.0509562 0.1665660 0.385187 0.573488 0.653394 0.676991
Mort 2 2016 1 2016_2_Mort 0.070000000 0.070000000 0.07000000 0.0700000 0.0700000 0.070000 0.070000 0.070000 0.070000
Keep 2 2016 1 2016_2_Keep 0.000127299 0.000619648 0.00294983 0.0133095 0.0514304 0.139040 0.239326 0.311748 0.365208
Dead 2 2016 1 2016_2_Dead 0.010470900 0.013230100 0.01803870 0.0306615 0.0694440 0.154574 0.251785 0.323324 0.377406
', header = TRUE, stringsAsFactors = FALSE, na.strings = "NA")
# EXPLOITATION table for 2016
# F_Method: 3 Continuous_F;_(NOTE:_F_report_adjusts_for_seasdur_but_each_fleet_F_is_annual)
# F_report_units: 1(F)/(Fmsy);_with_F=Exploit(bio)
# _ _ _ Bio Bio Num Bio
# _ _ _ 1 2 3 4
# Yr Seas F_report 1_TRAWL 2_FIX
# 2016 1 0.114781 0.0106242 0.00742963
Fvals <- c(0.01062420, 0.00742963)
my.Lsel <- Lsel_table[Lsel_table$Factor == "Lsel", !names(Lsel_table) %in% c("Yr", "Factor", "Sex", "Label")]
my.Dead <- Lsel_table[Lsel_table$Factor == "Dead", !names(Lsel_table) %in% c("Yr", "Factor", "Sex", "Label")]
my.Keep <- Lsel_table[Lsel_table$Factor == "Keep", !names(Lsel_table) %in% c("Yr", "Factor", "Sex", "Label")]
my.Mort <- Lsel_table[Lsel_table$Factor == "Ret", !names(Lsel_table) %in% c("Yr", "Factor", "Sex", "Label")]
my.Ret <- Lsel_table[Lsel_table$Factor == "Mort", !names(Lsel_table) %in% c("Yr", "Factor", "Sex", "Label")]
natural_mortality <- 0.25
# calculate survival probability after accounting for natural mortality and
# commercial mortality from both fleets using two approaches
# first approach
dead.beta <- Fvals * my.Dead[, 2:ncol(my.Dead)]
surv.prob1 <- exp(-(natural_mortality + colSums(dead.beta[,c('X49', 'X51', 'X53', 'X55', 'X57', 'X59', 'X61', 'X63', 'X65')])))
surv.prob1
# second approach
# commercial mortality * selectivity * (retain + (1 - retain) * discard)
comm.mort.beta <- Fvals * my.Lsel[, 2:ncol(my.Lsel)] * (my.Ret[, 2:ncol(my.Ret)] + (1 - my.Ret[, 2:ncol(my.Ret)]) * my.Mort[, 2:ncol(my.Mort)])
surv.prob2 <- exp(-(natural_mortality + colSums(comm.mort.beta[,c('X49', 'X51', 'X53', 'X55', 'X57', 'X59', 'X61', 'X63', 'X65')])))
surv.prob2
round(surv.prob1, 7) == round(surv.prob2, 7)
