Getting mortality from atlantisom: backing out F and M from Z

Introduction

This page documents initial testing of the atlantisom package in development at https://github.com/r4atlantis/atlantisom using three different Atlantis output datasets. Development of atlantisom began at the 2015 Atlantis Summit in Honolulu, Hawaii, USA. On this page we demonstrate use of atlantisom on both California Current (CCA) and Norwegian Barents Sea Atlantis (NOBA) output files to test function new functions that use calc_Z calculating total mortality to back out F and M. These are used initially for comparison with stock assessment model fishing mortality outputs, and to parameterize stock assessment natural mortality.

Previous comparisons of estimated total mortality with the output[scenario.name]Mort.txt file demonstrated what the Atlantis manual already tels us: >This file is currently only useful for looking at relative M vs F values for a species, as itdoes not give accurate mortalities

We can assumed that the scale of Z from this file will not be the same as our calculation, but since both F and M are rescaled, even from input F values, we need to back out the realized F to subtract from this Z to get M. This approach seems simpler than trying to estimate M directly using predation.

We will first use CCA here because the contrasting F scenario has been implemented. We will test with sardines only first.

All model setup and configuration is described here.

library(tidyr)
require(dplyr)
library(ggplot2)
library(data.table)
library(here)
library(ggforce)
library(ggthemes)
library(atlantisom)

initCCA <- TRUE
initNOBA <- FALSE

species_ss <- c("Pacific_sardine")

#function to make a config file? need one for each atlantis run

if(initCCA) source(here("config/CC3Config.R"))

if(initNOBA) source(here("config/NOBAaaConfig.R"))

#Load functional groups
funct.groups <- load_fgs(dir=d.name,
                         file_fgs = functional.groups.file)
#Get just the names of active functional groups
funct.group.names <- funct.groups %>% 
  filter(IsTurnedOn == 1) %>%
  select(Name) %>%
  .$Name

# default run_truth setup will save the file, so check for that first

if(!file.exists(file.path(d.name, 
                          paste0("output", scenario.name, "run_truth.RData")))){
  #Store all loaded results into an R object
  truth <- run_truth(scenario = scenario.name,
                     dir = d.name,
                     file_fgs = functional.groups.file,
                     file_bgm = box.file,
                     select_groups = funct.group.names,
                     file_init = initial.conditions.file,
                     file_biolprm = biol.prm.file,
                     file_runprm = run.prm.file
  )
} else{
  truth <- get(load(file.path(d.name,
                              paste0("output", scenario.name, "run_truth.RData"))))
}

Estimate total mortality with calc Z function

Trying to wrap these steps into something simpler.

# make a function for this
source(here("config/census_spec.R"))

Output timestep toutinc for the population is 73,

so steps per year in run_truth output is 5

and the number of output steps in run_truth output is 600.

Calculating Z and F at different levels of aggregation

Here we test outputs of Z for several levels of aggregation using files generated previously:

# add YOY file to the config files
YOY <- load_yoy(d.name, paste0("output", scenario.name, "YOY.txt"))

# load biolprm in some initialize file?
biol <- load_biolprm(d.name, biol.prm.file)

## Warning in `[<-.data.frame`(`*tmp*`, , 3:12, value = structure(list(`1` =
## structure(c(469L, : provided 13 variables to replace 10 variables

# get code matching species name to split YOY file
code_ss <- funct.groups$Code[which(funct.groups$Name == species_ss)]

# cut to a single species in YOY file
YOY_ss <- YOY %>%
  select(Time, paste0(code_ss, ".0"))

# numbers at agecl at full resolution (all polygons and layers)
truenums_ss <- truth$nums[truth$nums$species == species_ss,]

survey_testNall <- readRDS(file.path(d.name, paste0(scenario.name, "survey_testNall.rds")))

# numbers at agecl aggregated over layers, retain polygons (output of create_survey)
survey_testN_ss <- survey_testNall[survey_testNall$species == species_ss,] 

# numbers at agecl aggregated over layers and polygons
Natage_ss <- readRDS(file.path(d.name, paste0(scenario.name, "natage_census_sard_hake.rds")))

Natage_ss <- Natage_ss[Natage_ss$species == species_ss,]

#calc_Z <- function(yoy, nums, fgs, biolprm, toutinc)

fullresZ <- calc_Z(yoy = YOY_ss,
                   nums = truenums_ss,
                   fgs = funct.groups,
                   biolprm = biol,
                   toutinc = runpar$toutinc)

surveyresZ <- calc_Z(yoy = YOY_ss,
                   nums = survey_testN_ss,
                   fgs = funct.groups,
                   biolprm = biol,
                   toutinc = runpar$toutinc)

sampleresZ <- calc_Z(yoy = YOY_ss,
                   nums = Natage_ss,
                   fgs = funct.groups,
                   biolprm = biol,
                   toutinc = runpar$toutinc)

What is F at each timestep?

A relative annual F is found in outputMort.txt! Sadly, it is wrong.

file.mort <- file.path(d.name, paste0("output", scenario.name, "Mort.txt"))

mortish <- read.table(file.mort, header = TRUE)

relF_ss <- mortish %>%
  select(Time, relF = paste0(code_ss, ".F"))

relM_ss <- mortish %>%
  select(Time, relM = paste0(code_ss, ".M"))

Zish <- merge(relF_ss,relM_ss) %>%
  mutate(relZ = relF + relM)

plotF <- ggplot(data=relF_ss, aes(x=Time/365, y=relF)) +
  geom_line() +
  theme_tufte() +
  labs(subtitle = paste(scenario.name, species_ss))

plotF

Estimate F from catch and biomass at different levels of aggregation

Ideally we would want this calculation to be in numbers of individuals to calculate F. However we do not have reliable catch in numbers output from catch.nc for legacy models. Therefore I will build a switch in to check for legacy models to do this calculation in biomass, otherwise the default calculation should use catch.nc catch in numbers compared to the numbers (in the previous ouput timestep) in nums from truth.

For CCA, Catch in biomass at each catch output timestep is read in from catch.txt:

truecatchbio <- load_catch(d.name, file_catch = catch.file, fgs = funct.groups)

truecatchbio_ss <- truecatchbio[truecatchbio$species == species_ss,]

rm(truecatchbio)

Total biomass at each output timestep, or should this be average biomass? Or initial biomass? Isaac used initial biomass for the year from the biomass index and total catch for the year so comparing catch at day 365 to biomass at day 0 for that year.

Since we have only annual B then we can just read in true B from the biomass index:

atBtxt2 <- read.table(file.path(d.name, paste0("output", scenario.name, "BiomIndx.txt")), header=T)
groupslookup <- funct.groups %>%
  filter(IsTurnedOn > 0)

atBtxt2tidy <- atBtxt2 %>%
  select(Time:DIN) %>%
  #select(Time, FPL:DIN) %>%
  rename_(.dots=with(groupslookup, setNames(as.list(as.character(Code)), Name))) %>%
  gather(species, biomass, -Time) %>%
  filter(species %in% species_ss)

## Warning: rename_() is deprecated. 
## Please use rename() instead
## 
## The 'programming' vignette or the tidyeval book can help you
## to program with rename() : https://tidyeval.tidyverse.org
## This warning is displayed once per session.

Output time units match in these 2 files (days). I want to recode the time units so that year 1 is time 0 for biomass and time 365 for catch, then get annbioF:

annbio <- atBtxt2tidy %>%
  mutate(yr = (Time/365)+1)

anncatch <- truecatchbio_ss %>%
  mutate(yr = time/365)

annbioF <- merge(annbio, anncatch) %>%
  mutate(annbioF = atoutput/biomass) %>%
  arrange(yr)

This annual biomass based estimate of F is scaled down from the input F (show in points in the plot).

plotannbioF <- ggplot(data=annbioF, aes(x=yr, y=annbioF)) +
  geom_line() +
  theme_tufte() +
  labs(subtitle = paste(scenario.name, species_ss))

plotannbioF + geom_point(data=relF_ss, aes(x=Time/365, y=relF))

How to turn timestep Z into annual Z for proper comparison?

Two approaches: annual nums each year, calculate survival between years (minus annual recruitment), turn into annual Z, or add timestep Zs to get annual?

First an annual estimate based on end year numbers snapshots:

totnums <- aggregate(atoutput ~ species + time, data = truenums_ss, sum) %>%
    mutate(time.days = (time+1)*runpar$toutinc) %>% #makes time 0 into days 0->73, etc
    mutate(yr = ceiling(time.days/365))  # yr 1 is 0:stepsperyr to match recruits yr1
    #mutate(yr = floor(time.days/365)) # includes 0 value, yr 1 is 

#only want numbers at end of year, snapshot
totnumsann <- totnums %>%
  filter(time.days %in% seq(365, max(time.days), by=365))

# mg carbon converted to wet weight in tonnes
k_wetdry <- biol$kgw2d / 1000000000

# WARNING only works for CCA because YOY.txt rows 2:end are already in numbers
# and we are merging out the incorrect and irrelevant YOY row 1 (Time=0)
# also hardcoded for sardine example
recnums <- YOY_ss %>%
  mutate(yr = as.integer(round(YOY_ss$Time)/365)) %>%
  mutate(recnums = SAR.0/k_wetdry/biol$redfieldcn)

annZ1 <- merge(totnumsann, recnums) %>%
  mutate(recnums = replace_na(recnums, 0)) %>%
  mutate(numslessrec = atoutput - recnums) %>%
  mutate(surv = numslessrec/lag(atoutput, default = first(numslessrec))) %>%
  mutate(Z = -1 * log(surv))

Second summing within year Z from calc_Z:

annZ2 <- fullresZ %>%
  mutate(yr = floor(time/stepperyr)+1) %>%
  group_by(species, yr) %>%
  summarise(Z = sum(atoutput))

Compare methods:

plotZ <- ggplot() +
  geom_line(data=Zish, aes(x=Time/365, y=relZ, color="mort.txt Z")) +
  geom_point(data=annZ1, aes(x=Time/365, y=Z, color="endyr Z")) +
  geom_point(data=annZ2, aes(x=yr, y=Z, color="sum timestep Z")) +
  theme_tufte() +
  labs(subtitle = paste(scenario.name, species_ss))

plotZ + ylim(0, 2.0)

Both seem to be tracking each other reasonably, and also general trends from mort.txt output. Discussion with Beth suggests that mort.txt output is an initial calculation combining deaths due to predation with deaths due to fishing, but that there is rescaling of predation (and possibly fishing?) afterward this output stage such that mort.txt does not provide final mortality. If the rescaling varies by timestep, then matching general patterns rather than the full internannual variability is probably the best we can do, since we don’t know the full interannual variability.

We can probably conclude that the calc_Z function is working as correctly as possible now, given that it is an approximation no matter what, and the differences we see here are likley from temporal resolution of the Z estimate.

Now do the same thing with F

F based on numbers would have this option, but we cant for CC because we only have annual catch in biomass output. Subtract only from annual (endyr) Z.

Now get M from Z-F

So what is M? This is an approximate annual M as Z-F from the annual estimate of F derived above.

annbioF <- select(annbioF, yr, annbioF)

annM1 <- merge(annZ1, annbioF, all.x = T) %>%
  mutate(M = Z-annbioF)

annZ2_Time <- annZ2 %>%
  mutate(Time = yr*365)

annM2 <- merge(annZ2_Time, annbioF) %>%
  mutate(M = Z-annbioF)

plotM <- ggplot() +
  geom_line(data=Zish, aes(x=Time/365, y=relZ, color="mort.txt Z")) +
  geom_line(data=Zish, aes(x=Time/365, y=relM, color="mort.txt M")) +
  geom_point(data=annM1, aes(x=Time/365, y=M, color="M = endyr Z-F")) +
  geom_point(data=annM2, aes(x=yr, y=M, color="M = sum timestep Z-F")) +
  theme_tufte() +
  theme(legend.position = "bottom", legend.box = "horizontal") +
  scale_color_discrete(NULL) +
  guides(color = guide_legend(nrow = 1)) +
  labs(subtitle = paste(scenario.name, species_ss))

plotM + ylim(0, 2.0)

And our estimated M time series:

plotestM <- ggplot(data=annM2, aes(x=yr, y=M)) +
  geom_line() +
  theme_tufte() +
  labs(subtitle = paste(scenario.name, species_ss))

plotestM + ylim(0, 0.5)

We still have a mismatch between Z in numbers and F calculated based on biomass, with some negative M values resulting from this calculation. It has improved somewhat.

Is there a mismatch between biomass based F and numbers based Z? We can’t do this comparison for CCA because our catch numbers from the .nc file have the pre-December 2015 bug.

Also, this may indicate that we need to redo calc_Z as cohort-specific after all.

Test calc_Z NOBA

Let’s try with NOBA and compare biomss and numbers based F directly. I’ll use NOBA herring as a species similar to sardines.

initCCA <- FALSE
initNOBA <- TRUE

if(initNOBA) source(here("config/NOBAaaConfig.R"))

species_ss <- "Norwegian_ssh"

#Load functional groups
funct.groups <- load_fgs(dir=d.name,
                         file_fgs = functional.groups.file)
#Get just the names of active functional groups
funct.group.names <- funct.groups %>% 
  filter(IsTurnedOn == 1) %>%
  select(Name) %>%
  .$Name

#Get true NOBAaa
if(!file.exists(file.path(d.name, 
                          paste0("output", scenario.name, "run_truth.RData")))){
  #Store all loaded results into an R object
  truth <- run_truth(scenario = scenario.name,
                     dir = d.name,
                     file_fgs = functional.groups.file,
                     file_bgm = box.file,
                     select_groups = funct.group.names,
                     file_init = initial.conditions.file,
                     file_biolprm = biol.prm.file,
                     file_runprm = run.prm.file
  )
} else{
  truth <- get(load(file.path(d.name,
                              paste0("output", scenario.name, "run_truth.RData"))))
}

source(here("config/census_spec.R"))

Lets try just the Z part for NOBA:

# make a function for this
# add YOY file to the config files
YOY <- load_yoy(d.name, paste0("output", scenario.name, "YOY.txt"))

# load biolprm in some initialize file?
biol <- load_biolprm(d.name, biol.prm.file)

## Warning in `[<-.data.frame`(`*tmp*`, , 3:12, value = structure(list(`1` =
## structure(c(124L, : provided 20 variables to replace 10 variables

# get code matching species name to split YOY file
code_ss <- funct.groups$Code[which(funct.groups$Name == species_ss)]

# cut to a single species in YOY file
YOY_ss <- YOY %>%
  select(Time, paste0(code_ss, ".0"))

# numbers at agecl at full resolution (all polygons and layers)
truenums_ss <- truth$nums[truth$nums$species == species_ss,]

#calc_Z <- function(yoy, nums, fgs, biolprm, toutinc)

#oops, need to generalize calc_Z for subannual timesteps in YOY!
#or input YOY only at 0, 365, etc since the numbers repeat for timesteps 
YOY_ss <- YOY_ss %>%
  filter(Time %in% seq(0, max(Time), by=365))

fullresZ <- calc_Z(yoy = YOY_ss,
                   nums = truenums_ss,
                   fgs = funct.groups,
                   biolprm = biol,
                   toutinc = runpar$toutinc)

# compare as above with mort.txt output
file.mort <- file.path(d.name, paste0("output", scenario.name, "Mort.txt"))

mortish <- read.table(file.mort, header = TRUE)

relF_ss <- mortish %>%
  select(Time, relF = paste0(code_ss, ".F"))

relM_ss <- mortish %>%
  select(Time, relM = paste0(code_ss, ".M"))

Zish <- merge(relF_ss,relM_ss) %>%
  mutate(relZ = relF + relM)

totnums <- aggregate(atoutput ~ species + time, data = truenums_ss, sum) %>%
    mutate(time.days = (time+1)*runpar$toutinc) %>% #makes time 0 into days 0->73, etc
    mutate(yr = ceiling(time.days/365))  # yr 1 is 0:stepsperyr to match recruits yr1
    #mutate(yr = floor(time.days/365)) # includes 0 value, yr 1 is 

#only want numbers at end of year, snapshot
totnumsann <- totnums %>%
  filter(time.days %in% seq(365, max(time.days), by=365))

# mg carbon converted to wet weight in tonnes
k_wetdry <- biol$kgw2d / 1000000000

# WARNING only works for CCA because YOY.txt rows 2:end are already in numbers
# and we are merging out the incorrect and irrelevant YOY row 1 (Time=0)
# also hardcoded for sardine example
recnums <- YOY_ss %>%
  mutate(yr = as.integer(round(YOY_ss$Time)/365)) %>%
  mutate(recnums = SSH.0/k_wetdry/biol$redfieldcn)

annZ1 <- merge(totnumsann, recnums) %>%
  mutate(recnums = replace_na(recnums, 0)) %>%
  mutate(numslessrec = atoutput - recnums) %>%
  mutate(surv = numslessrec/lag(atoutput, default = first(numslessrec))) %>%
  mutate(Z = -1 * log(surv))

## Warning in log(surv): NaNs produced

annZ2 <- fullresZ %>%
  mutate(yr = floor(time/stepperyr)+1) %>%
  group_by(species, yr) %>%
  summarise(Z = sum(atoutput))

plotZ <- ggplot() +
  geom_line(data=Zish, aes(x=Time/365, y=relZ, color="mort.txt Z")) +
  #geom_point(data=fullresZ, aes(x=time/stepperyr, y=atoutput, color="each timestep Z")) +
  geom_point(data=annZ1, aes(x=yr-1, y=Z, color="annual Z")) +
  geom_point(data=annZ2, aes(x=yr-1, y=Z, color="sum timestep Z")) +
  theme_tufte() +
  labs(subtitle = paste(scenario.name, species_ss))

plotZ + ylim(-1, 10.0)

## Warning: Removed 13 rows containing missing values (geom_point).

This Z is a bad match. But in a different way from cod where Z was huge; here the mort.txt Z is smaller than the estimated Z. This makes me question the Z calculation again.

So does Atlantis scale M up from the input as well?

Or maybe spring spawning herring leave the model area and I’m not accounting for that (is this a run with broken migration?). No evidence that these species migrate in documentation.

Our estimate of Z is still from the actual numbers alive at the end of each timestep and the known input recruitment prior to any mortality, so I think it really is the best estimate we have.

Two methods for F; comparison

NOBA F from mort.txt

plotF <- ggplot(data=relF_ss, aes(x=Time/365, y=relF)) +
  geom_line() +
  theme_tufte() +
  labs(subtitle = paste(scenario.name, species_ss))

plotF

NOBA biomass based F

# catch from txt output
truecatchbio <- load_catch(d.name, file_catch = catch.file, fgs = funct.groups)

# species in this output is the code not the name, putting in the name
truecatchbio_ss <- truecatchbio %>%
  left_join(funct.groups, by = c("species"="Code")) %>% 
  mutate(species = ifelse(!is.na(Name), Name, Code)) %>%
  select(species, agecl, polygon, layer, time, atoutput) %>%
  filter(species == species_ss) %>%
  arrange(time)

## Warning: Column `species`/`Code` joining factor and character vector, coercing
## into character vector

rm(truecatchbio)

#biomass from txt output
atBtxt2 <- read.table(file.path(d.name, paste0("output", scenario.name, "BiomIndx.txt")), header=T)
groupslookup <- funct.groups %>%
  filter(IsTurnedOn > 0)

atBtxt2tidy <- atBtxt2 %>%
  select(Time:DIN) %>%
  #select(Time, FPL:DIN) %>%
  rename_(.dots=with(groupslookup, setNames(as.list(as.character(Code)), Name))) %>%
  gather(species, biomass, -Time) %>%
  filter(species %in% species_ss)

#merge and estimate F
annbio <- atBtxt2tidy %>%
  mutate(yr = (Time/365)+1)

anncatch <- truecatchbio_ss %>%
  mutate(yr = time/365)

annbioF <- merge(annbio, anncatch) %>%
  mutate(annbioF = atoutput/biomass) %>%
  arrange(yr)

This annual biomass based estimate of F is not scaled down from the input F (show in points in the plot), unlike in the CCA.

plotannbioF <- ggplot(data=annbioF, aes(x=yr, y=annbioF)) +
  geom_line() +
  theme_tufte() +
  labs(subtitle = paste(scenario.name, species_ss))

plotannbioF + geom_point(data=relF_ss, aes(x=Time/365, y=relF))

NOBA numbers based F

This should be numbers caught divided by total numbers (of fully selected fish) from the previous timestep.

truecatchnums <- create_fishery_subset(dat = truth$catch,
                                 time = timeall,
                                 species = species_ss,
                                 boxes = boxall)  

# apply default sample fish to get numbers at age aggregated over polygon
catch_numsss <- sample_fish(truecatchnums, effNhigh)

# save this, needed as input to calc_Z?
saveRDS(catch_numsss, file.path(d.name, paste0(scenario.name, "truecatch_nums_Nssh.rds")))

catch_numsss <- readRDS(file.path(d.name, paste0(scenario.name, "truecatch_nums_Nssh.rds")))

totnum <- readRDS(file.path(d.name, paste0(scenario.name, "Natageclcensus.rds"))) 

totnumss <- filter(totnum, species == species_ss)

# work in progress 
# align nums at agecl previous timestep with catch current timestep
# for each agecl: time-1 catch_numsss$atoutput / time totnumss$atoutput
catch_numsss <- catch_numsss %>%
  select(species, agecl, catchtime=time, catchnum=atoutput) %>%
  mutate(Ftime = catchtime - 1)

totnumss <- totnumss %>%
  select(species, agecl, time, popnum=atoutput) %>%
  mutate(Ftime = time)

#what about annual Z from totnumss catch curve?
# each timestep use lm(agecl~ log(atoutput))
library(dplyr)
library(broom) # to change model results into data frame
catchcurve1 <- totnumss %>% 
  group_by(time) %>%
  do(tidy(lm(log(popnum) ~ agecl, data = .))) %>% 
  filter(term == "agecl") 

# or do for the year averaging natage at each timestep?
catchcurve2 <- totnumss %>%
  mutate(yr = floor(time/stepperyr)+1) %>%
  group_by(yr, agecl) %>%
  summarise(avgpopnum = mean(popnum)) %>%
  do(tidy(lm(log(avgpopnum) ~ agecl, data = .))) %>% 
  filter(term == "agecl") 

altZ1 <- catchcurve1 %>%
  select(time, estimate) %>%
  mutate(yr = floor(time/stepperyr)+1) %>%
  group_by(yr) %>%
  summarise(Z = sum(-estimate))

altZ2 <- catchcurve2 %>%
  mutate(Z=-estimate) %>%
  select(yr, Z)

#another alternative tracking cohorts
# get rid of agecl 1 to ignore recruitment
# still need rec_time to know when they go to next agecl
# survival is 
#  within year 
# SSH actually have two ages per agecl so this wont work
# also, when I look at an age class numbers are INCREASING each timestep
# which makes no sense unless there is migration in?

#not done yet dont use this
timestepnumF <- merge(totnumss, catch_numsss) %>%
  mutate(stepnF = (catchnum/popnum)) %>%
  arrange(Ftime, agecl)

Try this for cod, herring don’t have much F

species_ss <- "North_atl_cod"

# get code matching species name to split YOY file
code_ss <- funct.groups$Code[which(funct.groups$Name == species_ss)]

# cut to a single species in YOY file
YOY_ss <- YOY %>%
  select(Time, paste0(code_ss, ".0"))

# numbers at agecl at full resolution (all polygons and layers)
truenums_ss <- truth$nums[truth$nums$species == species_ss,]

#calc_Z <- function(yoy, nums, fgs, biolprm, toutinc)

#oops, need to generalize calc_Z for subannual timesteps in YOY!
#or input YOY only at 0, 365, etc since the numbers repeat for timesteps 
YOY_ss <- YOY_ss %>%
  filter(Time %in% seq(0, max(Time), by=365))

fullresZ <- calc_Z(yoy = YOY_ss,
                   nums = truenums_ss,
                   fgs = funct.groups,
                   biolprm = biol,
                   toutinc = runpar$toutinc)

# compare as above with mort.txt output
file.mort <- file.path(d.name, paste0("output", scenario.name, "Mort.txt"))

mortish <- read.table(file.mort, header = TRUE)

relF_ss <- mortish %>%
  select(Time, relF = paste0(code_ss, ".F"))

relM_ss <- mortish %>%
  select(Time, relM = paste0(code_ss, ".M"))

Zish <- merge(relF_ss,relM_ss) %>%
  mutate(relZ = relF + relM)

annZ2 <- fullresZ %>%
  mutate(yr = floor(time/stepperyr)+1) %>%
  group_by(species, yr) %>%
  summarise(Z = sum(atoutput))

plotZ <- ggplot() +
  geom_line(data=Zish, aes(x=Time/365, y=relZ, color="mort.txt Z")) +
  geom_point(data=fullresZ, aes(x=time/stepperyr, y=atoutput, color="each timestep Z")) +
  geom_point(data=annZ2, aes(x=yr-1, y=Z, color="sum timestep Z")) +
  theme_tufte() +
  labs(subtitle = paste(scenario.name, species_ss))

plotZ + ylim(-1, 10.0)

## Warning: Removed 1 rows containing missing values (geom_point).

Now estimate cod F from biomass and compare:

# catch from txt output
truecatchbio <- load_catch(d.name, file_catch = catch.file, fgs = funct.groups)

# species in this output is the code not the name, putting in the name
truecatchbio_ss <- truecatchbio %>%
  left_join(funct.groups, by = c("species"="Code")) %>% 
  mutate(species = ifelse(!is.na(Name), Name, Code)) %>%
  select(species, agecl, polygon, layer, time, atoutput) %>%
  filter(species == species_ss) %>%
  arrange(time)

## Warning: Column `species`/`Code` joining factor and character vector, coercing
## into character vector

rm(truecatchbio)

#biomass from txt output
atBtxt2 <- read.table(file.path(d.name, paste0("output", scenario.name, "BiomIndx.txt")), header=T)
groupslookup <- funct.groups %>%
  filter(IsTurnedOn > 0)

atBtxt2tidy <- atBtxt2 %>%
  select(Time:DIN) %>%
  #select(Time, FPL:DIN) %>%
  rename_(.dots=with(groupslookup, setNames(as.list(as.character(Code)), Name))) %>%
  gather(species, biomass, -Time) %>%
  filter(species %in% species_ss)

#merge and estimate F
annbio <- atBtxt2tidy %>%
  mutate(yr = (Time/365)+1)

anncatch <- truecatchbio_ss %>%
  mutate(yr = time/365)

annbioF <- merge(annbio, anncatch) %>%
  mutate(annbioF = atoutput/biomass) %>%
  arrange(yr)

plotannbioF <- ggplot(data=annbioF, aes(x=yr, y=annbioF)) +
  geom_line() +
  theme_tufte() +
  labs(subtitle = paste(scenario.name, species_ss))

plotannbioF + geom_point(data=relF_ss, aes(x=Time/365, y=relF))

Cod F based on biomass isn’t that different from the F reported in mort.txt.

Cod M if we believe our Z estimate from calc_Z:

annZ2 <- fullresZ %>%
  mutate(yr = floor(time/stepperyr)+1) %>%
  group_by(species, yr) %>%
  summarise(Z = sum(atoutput))

annbioF <- select(annbioF, yr, annbioF)

# annM1 <- merge(annZ1, annbioF, all.x = T) %>%
#   mutate(M = Z-annbioF)

annZ2_Time <- annZ2 %>%
  mutate(Time = yr*365)

annM2 <- merge(annZ2_Time, annbioF) %>%
  mutate(M = Z-annbioF)

plotM <- ggplot() +
  geom_line(data=Zish, aes(x=Time/365, y=relZ, color="mort.txt Z")) +
  geom_line(data=Zish, aes(x=Time/365, y=relM, color="mort.txt M")) +
  #geom_point(data=annM1, aes(x=Time/365, y=M, color="M = endyr Z-F")) +
  geom_point(data=annM2, aes(x=yr, y=M, color="M = sum timestep Z-F")) +
  theme_tufte() +
  theme(legend.position = "bottom", legend.box = "horizontal") +
  scale_color_discrete(NULL) +
  guides(color = guide_legend(nrow = 1)) +
  labs(subtitle = paste(scenario.name, species_ss))

plotM + ylim(0, 10.0)