knitr::opts_chunk$set(echo = TRUE) require(here) require(kableExtra) require(dplyr) devtools::load_all("C:\\Users\\chris\\Documents\\GitHub\\jsonlite") devtools::load_all("C:\\Users\\chris\\Documents\\GitHub\\RMAS") devtools::install_github("cmlegault/ASAPplots", dependencies=FALSE) orig <- "C:/Users/chris/Documents/Github/age-structured-data/SS_simplelog/SS_simplelogistic/" out<- r4ss::SS_output(orig, covar=F, verbose=FALSE) MAS <- output_plots(data.dir = "data", years = 1987:2006, ages=seq(1,10), pop_name = "populations.txt", rep_name = "mas_report.txt", figs_dir="plots") in_control <- r4ss::SS_readctl_3.30(paste0(orig,"control.ss")) MAS_N<-matrix(as.numeric(unlist(strsplit(MAS$n_at_age," "))), byrow=T, ncol=20) SS_N <- out$natage[,c("Yr","Beg/Mid", "Era", min(out$agebins):(max(out$agebins)+2))] names(SS_N)[2] <- "pd" SS_inc <-SS_N %>% filter(pd=="B", Yr %in% seq(1987,2006)) %>% select(seq(4,12)) SS_plus <-SS_N %>% filter(pd=="B", Yr %in% seq(1987,2006)) %>% select(seq(13,15)) %>% rowSums() SS_N <- cbind(SS_inc, SS_plus) %>% t() r4ss::SS_plots(out, plot=c(1:6,8:24), html=FALSE, verbose = FALSE, uncertainty=FALSE)
This is a benchmark run with data generated using Stock Synthesis' bootstrap function on a species with life history parameters and weight-at-age similar to Pacific hake. However, fishery dynamics were greatly simplified. The model uses empirical weight-at-age for growth, a logistic age-based selectivity curve for the fishery and acoustic survey, and fixed values of natural mortality and steepness. Model years spanned r out$startyr
to r out$endyr
with r out$nseasons
season(s) and r out$nsexes
sex(es). Age bins were r out$agebins
. Other key parameters are given below in Table 1:
partable <- in_control$MG_parms[!grepl("F_|Recr|Impl|Early",rownames(in_control$MG_parms)),] kable(partable[,c("INIT","PHASE")], format="html", caption = "Table 1: Life History parameters used to generate data.") %>% kable_styling()
Both models kept most life history parameters fixed and estimated values for $R_0$, recruitment deviations, fishing, and selectivity parameters. Parameter estimates for SS and MAS are given below:
partable <- out$parameters %>% filter(!grepl("F_|Recr|Impl|Early",Label)) kable(partable[,c("Label","Value", "Phase")], format="html", caption = "Table 2: Life History parameters estimated by Stock Synthesis.") %>% kable_styling() parMAS <- MAS$parameters %>% filter(!grepl("fishing|recr",Name)) kable(parMAS, format="html", caption = "Table 2: Parameters estimates in MAS") %>% kable_styling()
We found a number of inconsistencies doing this comparison. First, to ensure SS and MAS values were comparable, we needed to implement specifying early recruitment deviations in MAS. Second, we made a number of changes to the way empirical weight-at-age was handled in MAS. Finally, we found the lognormal likelihood function specified in MAS was incorrect. After these were fixed, we were able to get reasonably consistent estimates for different indices.
knitr::include_graphics(paste0(orig,"plots/data_plot.png"))
knitr::include_graphics(paste0(orig,"plots/bio5_weightatsize.png"))
knitr::include_graphics("C:/Users/chris/Documents/GitHub/RMAS/plots/ObsVsExpectedSurvey IndexTotal.png") knitr::include_graphics("C:/Users/chris/Documents/GitHub/RMAS/plots/ObsVsExpectedCatch BiomassTotal.png")
rec_ss<- out$parameters %>% filter(grepl("RecrDev",Label)) %>% select(Value) rec_mas <- MAS$parameters %>% filter(grepl("recruitment",Name)) %>% select(Value) ind <- out$startyr:out$endyr plot(NA, xlim=c(out$startyr,out$endyr), ylim=c(-3,3), xlab="Year", ylab="Recruitment deviation") points(rec_ss[1:(length(ind)),"Value"] ~ ind) lines(as.numeric(as.character(rec_mas$Value))~ind)
F_ss<- out$parameters %>% filter(grepl("F_",Label)) %>% select(Value) F_mas <- MAS$parameters %>% filter(grepl("fishing",Name)) %>% select(Value) plot(NA, xlim=c(out$startyr,out$endyr), ylim=c(0,0.5), xlab="Year", ylab="F") points(F_ss[1:length(ind),"Value"] ~ ind) lines(as.numeric(as.character(F_mas$Value))~ind)
Points are estimated numbers from Stock Synthesis; lines are estimated numbers from MAS.
plot_NAtAge(SS_N, MAS_N, pdf_on=FALSE, byage = TRUE, years = 1987:2006, ages=seq(1,10))
Points are estimated numbers from Stock Synthesis; lines are estimated numbers from MAS.
plot_NAtAge(SS_N, MAS_N, pdf_on=FALSE, byage = F, years = 1987:2006, ages=seq(1,10))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.