Start with saved wrapper weight at age output for NOBA cod

source(here("config/NOBA2Config.R"))
NOBAom_cod <- readRDS(file.path(d.name, paste0(scenario.name, "omlist_ss.rds")))

fgs <- NOBAom_cod$funct.group_ss
biolprm <- NOBAom_cod$biol

wtage <- readRDS(file.path(d.name, paste0(scenario.name, "survObsWtAtAge.rds")))
wtage <- wtage[[1]] %>%
  filter(species == "North_atl_cod")

wageplot <- ggplot(wtage, aes(time, atoutput)) +
  geom_line(aes(colour = factor(agecl))) +
  theme_tufte() +
  theme(legend.position = "bottom") +
  xlab("model timestep (5 per year)") +
  ylab("average individual weight (g)") +
  labs(subtitle = paste0(scenario.name, " annual mid year sample"))

wageplot + facet_wrap(c("species"), scales="free_y")

If we were to do simple interpolation on this output, we can approximate weight at true age class:

# a function to go into atlantisom, sketched here, reuse calc_stage2age as possible

wtagecl <- wtage
annages <- NOBAom_cod$truenumsage_ss %>%
  filter(species == "North_atl_cod") %>%
  group_by(species, time, agecl) %>%
  rename(trueage = agecl) %>%
  summarise(truenatage = sum(atoutput))
  
#calc_wtTrueage <- function(wtagecl, annages, fgs) {
  # get number of ages per agecl from fgs.file
  #fgs <- atlantisom::load_fgs(dir=d.name, file_fgs = functional.groups.file)
  species.code <- fgs$Code
  #turnedon <- fgs[fgs$IsTurnedOn > 0, ]

  # Figure out the groups that have multiple ages classes in each stage (or
  # cohort) this is from the full model
  multiple_ages <- fgs[fgs$NumAgeClassSize>1, c(1,4,10)]
  
  # match species in wtagecl_data to those species in the fgs file
  sppwt <- levels(wtagecl$species)
  multiple_ages <- multiple_ages[multiple_ages$Name %in% sppwt, ]
  
  #now get indices
  names <- multiple_ages$Code
  num_multi_age <- dim(multiple_ages)[1]

  # timesteps from input wt file
  ntimesteps <- length(unique(wtagecl$time))

  # make a dataframe the right dimensions for all age classes
  # has species, time from wtage and trueage, truenatage from annages
  # add agecl column
  # finds the average age of each agecl each timestep based on truenatage
  wtage_out <- annages %>%
    dplyr::semi_join(wtagecl, by = c("species", "time")) %>%
    dplyr::left_join(multiple_ages, by = c("species" = "Name")) %>%
    dplyr::mutate(agecl = as.integer(ceiling(trueage/NumAgeClassSize))) %>%
    dplyr::group_by(species, time, agecl) %>%
    dplyr::mutate(avgage = weighted.mean(trueage, truenatage))
## Warning: Column `species` joining character vector and factor, coercing into
## character vector
  wtage_avgage <- wtage_out %>%
    select(species, time, agecl, avgage) %>%
    distinct()
  
  # find weight increment for a timestep
  # muweight for agecl+1 - muweight for agecl
  wtagecl_inc <- wtagecl %>%
    dplyr::arrange(species, time, agecl) %>%
    #dplyr::left_join(multiple_ages, by = c("species" = "Name")) %>%
    dplyr::group_by(species, time) %>%
    dplyr::mutate(increment = case_when(
      agecl==1 ~ atoutput,
      agecl>1 ~ atoutput - dplyr::lag(atoutput)
      )) %>%
    dplyr::left_join(wtage_avgage) %>%
    dplyr::mutate(avgageinc = case_when(
      agecl==1 ~ avgage,
      agecl>1 ~ avgage - dplyr::lag(avgage)
      )) 
## Joining, by = c("species", "agecl", "time")
## Warning: Column `species` joining factor and character vector, coercing into
## character vector
  # complex, but works for 2 ages/agecl, need to test more, need to fix oldest
  wtage_out <- wtage_out %>%
    dplyr::select(species, time, trueage, agecl, avgage) %>%
    dplyr::left_join(wtagecl_inc) %>%
    dplyr::group_by(species, time) %>%
    dplyr::mutate(wtIntage = case_when(
      trueage==1 ~ trueage/avgage*increment,
      (trueage>1 & trueage<avgage) ~ 
        (1-(avgage-trueage)/avgageinc)*increment + lag(atoutput),
      (trueage>1 & trueage>avgage) ~ 
        ((trueage-avgage)/lead(avgageinc))*lead(increment) + atoutput
        ))
## Joining, by = c("species", "time", "agecl", "avgage")
  # plot this--how bad would piecewise linear interpolation be?
  wtageclann <- ggplot(wtagecl_inc, aes(avgage, atoutput)) +
    geom_point() +
    geom_line(aes(colour = factor(time))) + 
    scale_x_continuous(minor_breaks = c(0:20)) +
    theme_tufte() +
    theme(legend.position = "none",
          panel.grid.minor.x = element_line(colour = "grey50"),
          panel.grid.major.x = element_line(colour = "grey50"))
  
  wtageclann + geom_point(data = wtage_out, 
             mapping = aes(trueage, wtIntage))
## Warning: Removed 144 rows containing missing values (geom_point).

  # maybe not so bad

  # fix oldest age

  #clean up to have only standard columns
  finaldata <- data.frame("species" = wtage_out$species,
                          "agecl" = wtage_out$trueage, "polygon" = NA, "layer" = NA,
                          "time" = wtage_out$time, "atoutput" = wtage_out$wtIntage)

  # output full set of muweight at newly created agecl

#}

Run wrappers and plot interpolated results for multiple NOBA species:

NOBAom <- om_init(here("config/NOBA2config.R"))

# 3 species: 40, 20, 6 age classes
NOBAom_redcodcap <- om_species(c("Redfish", "North_atl_cod", "Capelin"), NOBAom)

rm(NOBAom) #free some memory

Fix bug in omdimensions.R; need to vectorize:

omlist_ss <- NOBAom_redcodcap # multispecies call

n_age_classes <- omlist_ss$funct.group_ss$NumCohorts # works
age_classes <- 1:n_age_classes # nope, only first spp, needed in SS3 writing?

age_classes <- sapply(n_age_classes, seq)
names(age_classes)<-survspp

n_annages <- n_age_classes * omlist_ss$funct.group_ss$NumAgeClassSize # works
annages <- 1:n_annages # nope, only first spp, needed in SS3 writing?

annages <- sapply(n_annages, seq)
names(annages)<-survspp

# dimensioning survey in 3sppsurvey.R
survspp <- omlist_ss$species_ss # works
surveffic <- data.frame(species=survspp,
                     efficiency=rep(1.0,length(survspp)))  # works
survselex <- data.frame(species=rep(survspp, each=n_annages), # nope! 
                        agecl=rep(c(1:n_annages),length(survspp)),
                        selex=rep(1.0,length(survspp)*n_annages))
surveffN <- data.frame(species=survspp, effN=rep(1000, length(survspp))) # works
surv_cv <- data.frame(species=survspp, cv=rep(0.1,length(survspp))) # works

#dimensioned for annual ages rather than age classes
survselex <- data.frame(species=rep(survspp, n_annages), #  
                        agecl=unlist(sapply(n_annages,seq)),
                        selex=rep(1.0,sum(n_annages)))

New omdimensions file:

#survey species inherited from omlist_ss
survspp <- omlist_ss$species_ss
# survey season and other time dimensioning parameters
# generalized timesteps all models
noutsteps <- omlist_ss$runpar$tstop/omlist_ss$runpar$outputstep
timeall <- c(0:noutsteps)
stepperyr <- if(omlist_ss$runpar$outputstepunit=="days") 365/omlist_ss$runpar$toutinc
midptyr <- round(median(seq(0,stepperyr)))

# model areas, subset in surveyconfig
allboxes <- c(0:(omlist_ss$boxpars$nbox - 1))

# fishery output: learned the hard way this can be different from ecosystem outputs
fstepperyr <- if(omlist_ss$runpar$outputstepunit=="days") 365/omlist_ss$runpar$toutfinc

# survey selectivity (agecl based)
# should return all age classes fully sampled (Atlantis output is 10 age groups per spp)
n_age_classes <- omlist_ss$funct.group_ss$NumCohorts
# changed below for multiple species
age_classes <- sapply(n_age_classes, seq)
names(age_classes)<-survspp

n_annages <- n_age_classes * omlist_ss$funct.group_ss$NumAgeClassSize
# changed below for multiple species
annages <- sapply(n_annages, seq)
names(annages)<-survspp

New config files needed for multispecies pulls:

3sppsurvey.R configures the fishery independent survey and looks like this:

# Default survey configuration here is a census
# Need to define survey season, area, efficiency, selectivity

#Atlantis model timestep corresponding to the true output--now from census_spec.R
timestep <- stepperyr #5

#Which atlantis timestep does the survey run in?--now from census_spec.R
# with 5 output steps per year, 0 is Jan-Feb-midMar, 1 is midMar-Apr-May, 
# 2 is June-July-midAug, 3 is midAug-Sept-Oct, 4 is Nov-Dec (ish)

survey_sample_time <- midptyr #2; defined in omdimensions.R

#The last timestep to sample
total_sample <- noutsteps-1 #495

#Vector of indices of survey times to pull
survey_sample_full <- seq(survey_sample_time,
                          total_sample, by=timestep)

survtime <- survey_sample_full

# survey area
# should return all model areas
survboxes <- allboxes

# survey efficiency (q)
# should return a perfectly efficient survey 
surveffic <- data.frame(species=survspp,
                     efficiency=rep(1.0,length(survspp)))

# survey selectivity (agecl based)
# this is by age class, need to change to use with ANNAGEBIO output
#survselex <- data.frame(species=rep(survspp, each=n_age_classes),
#                     agecl=rep(c(1:n_age_classes),length(survspp)),
#                     selex=rep(1.0,length(survspp)*n_age_classes))

# for annage output
survselex <- data.frame(species=rep(survspp, n_annages), #  
                        agecl=unlist(sapply(n_annages,seq)),
                        selex=rep(1.0,sum(n_annages)))

# effective sample size needed for sample_fish
# this effective N is high but not equal to total for numerous groups
surveffN <- data.frame(species=survspp, effN=rep(1000, length(survspp)))

# survey index cv needed for sample_survey_xxx
# cv = 0.1
surv_cv <- data.frame(species=survspp, cv=rep(0.1,length(survspp)))

# length at age cv for input into calc_age2length function
# function designed to take one cv for all species, need to change to pass it a vector
lenage_cv <- 0.1

# max size bin for length estimation, function defaults to 150 cm if not supplied
maxbin <- 200

3sppfishery.R configures the fishery dependent data and looks like this:

# Default fishery configuration here is a census
# Inherits species from input omlist_ss
fishspp <- omlist_ss$species_ss 

#Number of years of data to pull
nyears <- 50

#Atlantis initialization period in years
burnin <- 30

# fishery output: learned the hard way this can be different from ecosystem outputs
fstepperyr <- if(omlist_ss$runpar$outputstepunit=="days") 365/omlist_ss$runpar$toutfinc


# same time dimensioning parameters as in surveycensus.R
#Vector of indices of catch in numbers to pull (by timestep to sum)
fish_sample_full <- c(0:total_sample)  #total_sample defined in sardinesurvey.R
fish_burnin <- burnin*fstepperyr+1
fish_nyears <- nyears*fstepperyr
fish_times <- fish_sample_full[fish_burnin:(fish_burnin+fish_nyears-1)]
fish_timesteps <- seq(fish_times[fstepperyr], max(fish_times), by=fstepperyr) #last timestep
#fish_years <- unique(floor(fish_times/fstepperyr)+1) # my original
fish_years <- unique(floor(fish_times/fstepperyr)) #from Christine's new sardine_config.R

fishtime <- fish_times


# fishery sampling area
# should return all model areas, this assumes you see everything that it caught
fishboxes <- c(0:(omlist_ss$boxpars$nbox - 1))

# effective sample size needed for sample_fish
# this effective N is high but not equal to total for numerous groups
fisheffN <- data.frame(species=survspp, effN=rep(1000, length(survspp)))

#this adjusts for subannual fishery output so original effN is for whole year
fisheffN$effN <- fisheffN$effN/fstepperyr 

# fishery catch cv can be used in sample_survey_biomass
# perfect observation
fish_cv <- data.frame(species=survspp, cv=rep(0.01,length(survspp)))

Fix multispecies bug in om_index function:

# original  
wtage <- data.frame(species=rep(survspp, each=max(age_classes)), #nope
                      agecl=rep(c(1:max(age_classes)),length(survspp)), #nope
                      wtAtAge=rep(1000.0,length(survspp)*max(age_classes))) #nope

# age_classes is now a list 1:n_age_classes by species
# this is still dimensioned by agecl not true age
wtage <- data.frame(species=rep(survspp, n_age_classes), 
                      agecl=unlist(sapply(n_age_classes,seq)), 
                      wtAtAge=rep(1000.0,sum(n_age_classes)))

Try index step again:

NOBAom_redcodcap_ind <- om_index(usersurvey = here("config/3sppsurvey.R"), 
                           userfishery = here("config/3sppfishery.R"),
                           omlist_ss = NOBAom_redcodcap, 
                           n_reps = 1, 
                           save = TRUE)

om_index works, om_comps runs for survey and produces age, length, weight output, produces fishery age output, but fails calling the age2length function.

check fishery observed age comp:

catch_age_comp <- readRDS(file.path(d.name, paste0(scenario.name, "fishObsAgeComp.rds")))
catch_age_comp <- catch_age_comp[[1]]

filter(catch_age_comp, species=="North_atl_cod")
filter(catch_age_comp, species=="Redfish")
filter(catch_age_comp, species=="Capelin")

#snippets from calc_age2length
nums <- filter(catch_age_comp, species=="Capelin")
#calculate mean length at age
mulen <- nums
muweight <- nums

#needed stuff from om_comps
#Get catch weights for length comp calc
  # aggregate true resn per fishery subset design
  catch_aggresnss <- atlantisom::aggregateDensityData(dat = omlist_ss$trueresn_ss,
                                                      time = fishtime,
                                                      species = survspp,
                                                      boxes = fishboxes)

  # aggregate true structn fishery subsetdesign
  catch_aggstructnss <- atlantisom::aggregateDensityData(dat = omlist_ss$truestructn_ss,
                                                         time = fishtime,
                                                         species = survspp,
                                                         boxes = fishboxes)

  #dont sample these, just aggregate them using median
  structn <- atlantisom::sample_fish(catch_aggstructnss, fisheffN, sample = FALSE)

  resn <- atlantisom::sample_fish(catch_aggresnss, fisheffN, sample = FALSE)
  
  # more snippets from calc_age2length
  # extract rows from structn and resn that match the rows in nums, which are only non-zeroes
resn.ind <- with(resn,paste(species,'.',agecl,'.',polygon,'.',layer,'.',time,sep=""))
num.ind <- with(nums,paste(species,'.',agecl,'.',polygon,'.',layer,'.',time,sep=""))
pick <- match(num.ind,resn.ind)
SRN <- resn$atoutput[pick] + structn$atoutput[pick]

# get weight-length parameters
li_a_use <- biolprm$wl[match(fgs$Code[match(nums$species,fgs$Name)],biolprm$wl[, 1]), 2]
li_b_use <- biolprm$wl[match(fgs$Code[match(nums$species,fgs$Name)],biolprm$wl[, 1]), 3]

#calc mean length and weight at age
mulen$atoutput <- ((biolprm$kgw2d*biolprm$redfieldcn*SRN)/(1000*li_a_use))^(1/li_b_use)
muweight$atoutput <- li_a_use*mulen$atoutput^li_b_use

upper.bins <- 1:maxbin
lower.bins <- c(0,upper.bins[-length(upper.bins)])
lenfreq <- NULL

# well, match does return NA for capelin, and this is what breaks the comparison using 
# the max function (returns NA). adding na.rm=TRUE to line 86
if(maxbin<max(mulen$atoutput)){
  print("Warning: maximum bin size is smaller than the longest fish in the sample. Fish above the maximum bin size will be removed from length compositions.")
}

groups <- unique(as.factor(structn$species)) 
CVlenage = 0.1
# but will NAs break the following?
for (irow in 1:nrow(mulen))
#for (irow in 1:500)
  {
  group <- nums$species[irow]
  igroup <- which(groups==group)
  box <- nums$polygon[irow]
  layer <- nums$layer[irow]
  time <- nums$time[irow]
  age <- nums$agecl[irow]

  sigma = sqrt(log((CVlenage^2)+1))
  muuse <- log(mulen$atoutput[irow]) - 0.5*(sigma^2)
  CumFracperbin <- plnorm(upper.bins,muuse,sigma)
  Fracperbin <- c(CumFracperbin[1],diff(CumFracperbin))
  natlength = Fracperbin*nums$atoutput[irow]
  results <- cbind(mulen[irow,],lower.bins,upper.bins, row.names = NULL) # add row.names = NULL?
  results$atoutput <- round(natlength,0)
  #results <- results[results$atoutput>0,]
  lenfreq <- rbind(lenfreq,results) # add row.names = NULL?
}

remove.zeroes = TRUE
#get rid of zero rows.
if (remove.zeroes) lenfreq <- lenfreq[lenfreq$atoutput>0,]

#output. natlength now has additional columns because of need to store length bin info
lenout <- NULL
lenout$mulen <- mulen
lenout$muweight <- muweight
lenout$natlength <- lenfreq

# ok, works but adds NA for Capelin age classes where no numbers were sampled.
# could crank up sample size or just live with this.

Try om_comps step again (now works by allowing NA in max comparison of calc_age2length():

NOBAom_redcodcap_comp <- om_comps(usersurvey = here("config/3sppsurvey.R"), 
                           userfishery = here("config/3sppfishery.R"),
                           omlist_ss = NOBAom_redcodcap, 
                           n_reps = 1, 
                           save = TRUE)

What do the interpolated weight at ages look like? Diagnostic plots from function look the same as above for cod (good), blank for Capelin (good? didnt need extrapolation but may have other problems), and bad for Redfish between true ages 2-15; lots of overestimation.

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

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

NOBAom_redcodcap <- readRDS(file.path(d.name, paste0(scenario.name, "omlist_ss.rds")))
# plot portion, just redfish, capelin is blank and cod looks the same
#library(ggforce)
survObsFullWtAtAge[[1]][[2]] + facet_wrap_paginate(~species, ncol=1, nrow = 1, page=3, scales = "free")

Compare ObsWtAtAge with ObsFullWtAtAge outputs by species:

# dataframe portion
intwtCod <- survObsFullWtAtAge[[1]][[1]] %>%
  filter(species == "North_atl_cod")

intwtRed <- survObsFullWtAtAge[[1]][[1]] %>%
  filter(species == "Redfish")

intwtCap <- survObsFullWtAtAge[[1]][[1]] %>%
  filter(species == "Capelin") #empty, no need for interpolation, but make intelligent

wtageCod <- survObsWtAtAge[[1]] %>%
  filter(species == "North_atl_cod")

wtageRed <- survObsWtAtAge[[1]] %>%
  filter(species == "Redfish")

wtageCap <- survObsWtAtAge[[1]] %>%
  filter(species == "Capelin")

# handmade plots:

wtagecl <- wtageCap

  wtagep <- ggplot(wtagecl, aes(agecl, atoutput)) +
    geom_point() +
    geom_line(aes(colour = factor(time))) + 
    scale_x_continuous(minor_breaks = c(0:max(wtagecl$agecl))) +
    theme_tufte() +
    theme(legend.position = "none",
          panel.grid.minor.x = element_line(colour = "grey50"),
          panel.grid.major.x = element_line(colour = "grey50"))

wtagep

Capelin is fine with WtAtAge output because it has age classes = true ages! Maybe just warn that a species with fgs$NumAgeClassSize = 1 will have NA output of interpolation functions, which will run if there is output for that species in the ANNAGEBIO.nc and ANNAGECATCH.nc files.

Redfish appear to have decreasing average weight at age as they get older:

wtagecl <- wtageRed

  wtagep <- ggplot(wtagecl, aes(agecl, atoutput)) +
    geom_point() +
    geom_line(aes(colour = factor(time))) + 
    scale_x_continuous(minor_breaks = c(0:max(wtagecl$agecl))) +
    theme_tufte() +
    theme(legend.position = "none",
          panel.grid.minor.x = element_line(colour = "grey50"),
          panel.grid.major.x = element_line(colour = "grey50"))

wtagep

Diagnose Redfish problem (it was not negative increments, but inadequate lagging and leading):

annageRed <- NOBAom_redcodcap$truenumsage_ss %>%
  filter(species == "Redfish")

annageCod <- NOBAom_redcodcap$truenumsage_ss %>%
  filter(species == "North_atl_cod")

annageCap <- NOBAom_redcodcap$truenumsage_ss %>%
  filter(species == "Capelin")

annages <- annageRed
wtagecl <- wtageRed
fgs <- NOBAom_redcodcap$funct.group_ss

  annages <- annages %>%
    group_by(species, time, agecl) %>%
    rename(trueage = agecl) %>%
    summarise(truenatage = sum(atoutput))

  # assuming fgs is data object already been read in by load_fgs elsewhere
  # Figure out the groups that have multiple ages classes in each stage (or
  # cohort)
  multiple_ages <- fgs[fgs$NumAgeClassSize>1, c(1,4,10)]

  # match species in wtagecl_data to those species in the fgs file
  sppwt <- levels(wtagecl$species)
  multiple_ages <- multiple_ages[multiple_ages$Name %in% sppwt, ]

  # make a dataframe the right dimensions for all age classes
  # has species, time from wtage and trueage, truenatage from annages
  # add agecl column
  # finds the average age of each agecl each timestep based on truenatage
  wtage_out <- annages %>%
    dplyr::semi_join(wtagecl, by = c("species", "time")) %>%
    dplyr::inner_join(multiple_ages, by = c("species" = "Name")) %>%
    dplyr::mutate(agecl = as.integer(ceiling(trueage/NumAgeClassSize))) %>%
    dplyr::group_by(species, time, agecl) %>%
    dplyr::mutate(avgage = weighted.mean(trueage, truenatage))
## Warning: Column `species` joining character vector and factor, coercing into
## character vector
  wtage_avgage <- wtage_out %>%
    select(species, time, agecl, avgage) %>%
    distinct()
  
  # find weight increment for a timestep
  # muweight for agecl+1 - muweight for agecl
  wtagecl_inc <- wtagecl %>%
    dplyr::arrange(species, time, agecl) %>%
    #dplyr::left_join(multiple_ages, by = c("species" = "Name")) %>%
    dplyr::group_by(species, time) %>%
    dplyr::mutate(increment = case_when(
      agecl==1 ~ atoutput,
      agecl>1 ~ atoutput - dplyr::lag(atoutput)
    )) %>%
    dplyr::left_join(wtage_avgage) %>%
    dplyr::mutate(avgageinc = case_when(
      agecl==1 ~ avgage,
      agecl>1 ~ avgage - dplyr::lag(avgage)
    ))
## Warning: Column `species` joining factor and character vector, coercing into
## character vector
  # complex, but works for 2 ages/agecl, need to test more, need to fix oldest
  # for more than 2 ages/agecl need to have a bigger lead/lag
  # achieved by faking a scalar for lead/lag n using max() function
  # should be ceiling I think for odd NumAgeClassSize, but not sure?
  # for oldest age(s), perhaps default=last(x) in lead statement
  wtage_out <- wtage_out %>%
    dplyr::select(species, time, trueage, NumAgeClassSize, agecl, avgage) %>%
    dplyr::left_join(wtagecl_inc) %>%
    dplyr::group_by(species, time) %>%
    dplyr::mutate(wtIntage = case_when(
      (agecl==1 & trueage<avgage) ~ (1-(avgage-trueage)/avgageinc)*increment,
      (agecl>1 & trueage<avgage) ~
        (1-(avgage-trueage)/avgageinc)*increment +
        dplyr::lag(atoutput, ceiling(max(NumAgeClassSize)/2)),
      (trueage>avgage) ~
        (trueage-avgage)/dplyr::lead(avgageinc, ceiling(max(NumAgeClassSize)/2), default=last(avgageinc))*dplyr::lead(increment, ceiling(max(NumAgeClassSize)/2), default=last(increment)) + atoutput
    ))


  # diagnostic plot--are interpolations where we expect them?
  wtageclann <- ggplot(wtagecl_inc, aes(avgage, atoutput)) +
    geom_point() +
    geom_line(aes(colour = factor(time))) +
    scale_x_continuous(minor_breaks = c(0:max(wtage_out$trueage))) +
    theme_tufte() +
    theme(legend.position = "none",
          panel.grid.minor.x = element_line(colour = "grey50"),
          panel.grid.major.x = element_line(colour = "grey50"))

  diag.p <- wtageclann + geom_point(data = wtage_out,
                          mapping = aes(trueage, wtIntage))

  # #clean up to have only standard columns
  # finaldata <- data.frame("species" = wtage_out$species,
  #                         "agecl" = wtage_out$trueage, "polygon" = NA, "layer" = NA,
  #                         "time" = wtage_out$time, "atoutput" = wtage_out$wtIntage)

  diag.p

Try the updated function again, any bugs with multiple species? Is it getting the weight at max age for both? OK so far.

omlist_ss <- NOBAom_redcodcap

interp_survWtAtAge <- calc_avgwtstage2age(wtagecl = survObsWtAtAge[[1]],
                                                     annages = omlist_ss$truenumsage_ss,
                                                     fgs = omlist_ss$funct.group_ss)

interp_survWtAtAge[[2]] + facet_wrap_paginate(~species, ncol=1, nrow = 1, page=1, scales = "free")

interp_survWtAtAge[[2]] + facet_wrap_paginate(~species, ncol=1, nrow = 1, page=2, scales = "free")

interp_survWtAtAge[[2]] + facet_wrap_paginate(~species, ncol=1, nrow = 1, page=3, scales = "free")

Extrapolating the oldest wt at age is now linear from the last agecl weight at age using default=last() in the dplyr::lead() calls.