APAM: American plaice assessment model

APAM is an age-based state-space stock assessment model developed for the American plaice stock on the Grand Bank of Newfoundland (NAFO divisions 3LNO). The base formulation is described in Perreault et al., 2020 and can be run by following the steps below. We note that this model formulation does not include data for the Spanish surveys due to confidentiality issues. How to fit some model variations (e.g. parameter map, M assumption) is also explained below; see the function documentation and Perreault et al., 2020 for additional details.

We also include details on how to run the M diagnostics from Perreault & Cadigan (2021, in review).

This is still a work in progress, and we aim to add other model diagnostics (e.g. retros, self sims), additional plots and tables and testing soon!

Installation

The latest version can be installed via:

#note this may produce warnings from RcppEigen, but does not affect model fitting/results
#see note at end of document on how to quiet warnings if interested 
devtools::install_github("SineAndie/APAM", dependencies=TRUE)

To run the base model

Running APAM consists of four steps:

1. reading in the data,
2. making the parameter map,
3. defining parameters and
4. model fitting.

Notes on how to modify inputs and model formulation are detailed in the next section.

library(APAM)

#read in data and structure for use in APAM
tmb.data<-make.tmb.data()

#make parameter map 
pmap<-make.map(tmb.data)

#prepare parameters
params<-make.parm(tmb.data, pmap)

#Finally, fit the model! (note: takes approx. 3 mins to run)
fit<-make.fit(tmb.data, pmap, params)

Outputs can be viewed with the make.plot function. There are a variety of plots available, and a few examples are shown below.

plots<-make.plots(fit)

#estimated population processes
plots$pop_process

#observed vs predicted survey indices ages 1-7
plots$index_fit1

#survey residuals
plots$resid_index

Explore model variations

The natural mortality assumption, including increasing M for years 1989-1996, can be changed within the make.tmb.data function. For example:

#turn off M split 
tmb.data_nosplit <- make.tmb.data(M.split = F) 

#change M assumption to 0.5 for all ages and years
M.matrix = matrix(0.5,nrow = tmb.data$Y, ncol = tmb.data$A)

tmb.data_newM <- make.tmb.data(M.matrix = M.matrix) 

Additionally, the process errors can be turned on/off through the make.map and make.parm functions.

pmap2<-make.map(tmb.data,no.pe = TRUE)
params2<-make.parm(tmb.data, pmap2, no.pe = TRUE)

#Then fit the model like usual
fit2<-make.fit(tmb.data, pmap2, params2)

#check for convergence/output nll 
fit$opt$message
#> [1] "relative convergence (4)"
fit$obj$fn()
#> [1] 3515.087
#> attr(,"logarithm")
#> [1] TRUE

fit2$opt$message
#> [1] "relative convergence (4)"
fit2$obj$fn()
#> [1] 3727.348
#> attr(,"logarithm")
#> [1] TRUE

The parameter map can be manually changed via setmap (still needs a bit of work/cleanup but functionality is there)

##default formulation: 
#setmap <- list(
    #   meanF= c("5","6"),
    #   stdF = c("5",rep("6+",length =   tmb.data$A-5)),
    #   ageFall = c("1",rep("2-11",10),rep("12-15",4)),
    #   ageSpring = c("1","2",rep("3-13",11),rep("14-15",2)),
    #   ageSpanish = NULL,
    #   stdcrl  = c(rep("5-6",2),rep("7-11",5),rep("12-14",3)),
    #   stdpe = rep("all", tmb.data$A-1),
    #   mapq = c(1:7,rep(NA,length = (tmb.data$A-1)-7))
    # )

#change to one meanF parameter and one std F parameter
setmap <- list(
  meanF= c("6"),
  stdF = c(rep("5+",length =   tmb.data$A-4)),
  ageFall = c("1",rep("2-11",10),rep("12-15",4)),
  ageSpring = c("1","2",rep("3-13",11),rep("14-15",2)),
  ageSpanish = NULL,
  stdcrl  = c(rep("5-6",2),rep("7-11",5),rep("12-14",3)),
  stdpe = rep("all", tmb.data$A-1),
  mapq = c(1:7,rep(NA,length = (tmb.data$A-1)-7))
)

pmap3<-make.map(tmb.data, setmap=setmap)

Natural mortality diagnostics

We can also run the natural mortality diagnostics from Perreault & Cadigan (2021, in review). To get profile likelihoods:

#to calculate profile likelihoods (takes approx 2 hrs)
profile <- make.profile(fit)

#to plot profile likelihood results 
prof_plots <- make.profile.plots(profile)

and local influence (LI) diagnostics. Note, in make.fit, make sure that do.hess=T or the LI diagnostics will return an error.

#to calculate age group LI slopes for the full model and data components (takes approx 15 mins)
LI_age <- make.LI(fit, age = T)

#plot the results
LI_age_plot <- make.LI.plots(LI_age)

#to calculate year group LI slopes (takes approx 45 mins)
LI_year <- make.LI(fit, year = T)

#plot the results
LI_year_plot <- make.LI.plots(LI_year)

#to get individual age and year LI slopes.
#note: for now this only returns the influence slopes from the full model fit (not data components since 
#the run time is quite long; takes approx 2 hrs to run full model)
LI <- make.LI(fit)

LI_plot <- make.LI.plots(LI)

Can also check the curvature to assess the suitability on the local influence diagnostics. Note, we show the code for the age groups below but similar for year group perturbations (takes approx 1.2 hrs) and individual perturbations (takes approx 14 hrs).

#check curvature (takes approx 20 mins to run)
curv_age <- make.curv(fit,LI_age, age = T)

#plot the results
curv_age_plot <- make.curv.plots(curv_age)

How to quiet RcppEigen warnings

Note: this is the workaround that has worked for me on Windows, but am open to any advice on more efficient/better methods.

1. find your R home directory via path.expand("~") (i.e. type this in the R console).
2. create a folder in the home directory called .R/; e.g. my home directory is C:/Users/myname/Documents, so I create C:/Users/myname/Documents/.R/
3. create a text file in the .R folder called Makevars.win
4. in Makevars.win writeCXXFLAGS = -Wno-ignored-attributes
5. save and install package - most warnings should be quiet now!