a4amanual.Rnw

\documentclass[a4paper,english,11pt]{article}
\usepackage{a4a}

\begin{document}

\title{Stock assessment with the \aFa statistical catch-at-age framework}

\input{authors}

\maketitle

\begin{abstract}
This chapter presents the statistical catch-at-age stock assessment model developed as part of the Assessment For All (\aFa) initiative of the European Commission Joint Research Centre. The stock assessment model framework is a non-linear catch-at-age model implemented in \href{http://www.r-project.org/}{\R}, \href{http://www.flr-project.org/}{\FLR} and \href{http://www.admb-project.org/}{\ADMB}, that can be applied rapidly to a wide range of situations with low parametrization requirements. The model structure is defined by submodels, which are the different parts that require structural assumptions. There are 5 submodels in operation: F-at-age, abundance indices catchability-at-age, recruitment, observation variance of catch-at-age and abundance indices-at-age, and population's age structurein the first year. All submodels use the same type of specification process, the \R formula interface, wich gives lot's of flexibility to explore models and combination of submodels.
\end{abstract}

\newpage
\tableofcontents
\newpage

<<knitr_opts, echo=FALSE, message=FALSE, warning=FALSE>>=
library(knitr)
library(formatR)
#thm = knit_theme$get("bclear") #moe, bclear
#knit_theme$set(thm)
opts_chunk$set(dev='png', cache=TRUE, fig.align='center', warning=FALSE, message=FALSE, dev.args=list(type="cairo"), dpi=96, highlight=TRUE, background='#F2F2F2', fig.lp="fig:", fig.pos="H", width=70, tidy=TRUE, out.width='.9\\linewidth')
# lattice theme
@

\section{Before starting}

\subsection{License, documentation and development status}

The software is released under the \href{https://joinup.ec.europa.eu/community/eupl/home}{EUPL 1.1}.

For more information on the \aFa methodologies refer to \href{http://icesjms.oxfordjournals.org/content/early/2014/04/03/icesjms.fsu050.abstract}{Jardim, et.al, 2014}, \href{http://icesjms.oxfordjournals.org/content/early/2014/03/31/icesjms.fsu043.abstract}{Millar, et.al, 2014} and \href{http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0154922}{Scott, et.al, 2016}.

Documentation can be found at \url{http://flr-project.org/FLa4a}. You are welcome to:

\begin{itemize}
	\item Submit suggestions and bug-reports at: \url{https://github.com/flr/FLa4a/issues}
	\item Send a pull request on: \url{https://github.com/flr/FLa4a/}
	\item Compose a friendly e-mail to the maintainer, see `packageDescription('FLa4a')`
\end{itemize}

\subsection{Installing and loading libraries}

To run the \code{FLa4a} methods the reader will need to install the package and its dependencies and load them. Some datasets are distributed with the package and as such need to be loaded too.

<<eval=FALSE>>=
# from CRAN
install.packages(c("copula","triangle", "coda", "grid", "gridExtra", "latticeExtra"))
# from FLR
install.packages(c("FLCore", "FLa4a"), repos="http://flr-project.org/R")
@

<<load_libraries, message=FALSE, warning=FALSE>>=
# libraries
library(devtools)
library(FLa4a)
library(XML)
library(reshape2)
library(ggplotFL)
# datasets
data(ple4)
data(ple4.indices)
data(ple4.index)
data(rfLen)
@

<<>>=
packageVersion("FLCore")
packageVersion("FLa4a")
@

\subsection{How to read this document}

The target audience for this document are readers with some experience in R and some background on stock assessment.

The document explains the approach being developed by \aFa for fish stock assessment and scientific advice. It presents a mixture of text and code, where the first explains the concepts behind the methods, while the last shows how these can be run with the software provided. Moreover, having the code allows the reader to copy/paste and replicate the analysis presented here.

The sections and subsections are as independent as possible, so it can be used as a reference document for the \code{FLa4a}. 

\subsection{Notation}

Along this chapter the notation presented in Table~\ref{tab:mathsnotation} will be used. Mathematical descriptions will be kept as simple as possible for readability.

\begin{table}
  \caption{Mathematical notation}
  \label{tab:mathsnotation}

  \begin{center}
    \begin{tabular}{cc}
      \hline
      Symbol & Description \\
      \hline
      \hline
       Variables & \\
       $C$ & catches \\
       $F$ & fishing mortality \\
       $M$ & natural mortality \\
       $R$ & recruitment \\
       $Q$ & vessel or fleet catchability \\
       $w$ & weights \\
       $l$ & likelihood \\
       $I$ & abundance index \\
       $S$ & spawning stock biomass \\
       $CV$ & coefficient of variation \\
       $D$ & residuals or deviances \\
       $N$ & normal distribution \\
       $\beta$ & parameter \\
       $a$ & stock-recruitment parameter \\
       $b$ & stock-recruitment parameter \\
       $\sigma^2$ & variance of catch \\
       $\tau^2$ & variance of index \\
       $\phi^2$ & variance of predicted recruitment \\
       $\upsilon^2$ & variance of residuals \\
       Subscripts & \\
       $a$ & age \\
       $y$ & year \\
       $C$ & catch \\
       $I$ & abundance index \\
       $N$ & normal distribution \\
       $s$ & survey \\
       $SR$ & stock recruitment relationship \\
       Superscripts and accents & \\
       $\hat{}$ & observation \\
       $\tilde{}$ & prediction \\
       $c$ & catches \\
       $s$ & abundance index \\
      \hline
    \end{tabular}
  \end{center}
\end{table}

\section{Introduction}

The \aFa stock assessment framework is based in a non-linear catch-at-age model implemented in \R, \FLR and \ADMB that can be applied rapidly to a wide range of situations with low setup requirements.

The framework is built of submodels which define the different parts of a statistical catch at age model that require structural assumptions. In the \aFa framework these are fishing mortality-at-age, abundance indicies catchability-at-age, recruitment, observation variances of catch-at-age and abundance indices-at-age, and abundance-at-age in the first year of the data series (see section~\ref{sec:math} for details).

Other important processes, like natural mortality, individual growth and reproduction, are treated as fixed(inputs??), as it's common in stock assessment methods. Nevertheless, the \aFa framework provides methods to condition these processes prior to the model fit, and propagate their uncertainty into the assessment process. See chapters XX and section XX.

The submodels formulation uses linear models, which opens the possibility of using the linear modelling tools available in \R. For example, \href{http://cran.r-project.org/web/packages/mgcv/index.html}{\code{mgcv}} gam formulas or factorial design formulas using \code{lm()}.

The 'language' of linear models has been developing within the statistical community for many years, and constitutes an elegant way of defining models without going through the complexity of mathematical representations. This approach makes it also easier to communicate among scientists:
\begin{itemize}
  \item \href{http://rspa.royalsocietypublishing.org/content/283/1393/147.short}{1965 J. A. Nelder}, notation for randomized block design
  \item \href{http://www.jstor.org/stable/info/2346786}{1973 Wilkinson and Rodgers}, symbolic description for factorial designs
  \item \href{http://books.google.com/books?isbn=0412343908}{1990 Hastie and Tibshirani}, introduced notation for smoothers
  \item \href{http://books.google.com/books?isbn=041283040X}{1991 Chambers and Hastie}, further developed for use in S
\end{itemize}

\section{Stock assessment framework maths description}\label{sec:math}

The stock assessment model behind the framework is based in two main equations to link observations to processes, the Baranov's catch equation and abundance indices equation.

Catches in numbers by age and year are defined in terms of the three quantities: natural mortality, fishing mortality and recruitment; using a modified form of the well known Baranov catch equation:

$$C_{ay} = \frac{F_{ay}}{F_{ay}+M_{ay}}\left(1 - e^{-(F_{ay}+M_{ay})}\right) R_{y}e^{-\sum (F_{ay} + M_{ay})} $$

Survey indices by age and year are defined in terms of the same three quantities with the addition of survey catchability:

$$I_{ays} = Q_{ays} R_{y}e^{-\sum (F_{ay} + M_{ay})}$$

Observed catches and observed survey indices are assumed to be log-normally distributed, or equivalently, normally distributed on the log-scale, with specific observation variance:

$$ \log \hat{C}_{ay} \sim \text{Normal} \Big( \log C_{ay}, \sigma^2_{ay}\Big) $$
$$ \log \hat{I}_{ays} \sim \text{Normal} \Big( \log I_{ays}, \tau^2_{ays} \Big) $$

The log-likelihood can now be defined as the sum of the log-likelihood of the observed catches:

$$  \ell_C = \sum_{ay} w^{(c)}_{ay}\ \ell_N \Big( \log C_{ay}, \sigma^2_{ay} ;\ \log \hat{C}_{ay} \Big) $$

and the log-likelihood of the observed survey indices as:

$$  \ell_I = \sum_s \sum_{ay} w^{(s)}_{ays}\ \ell_N \Big( \log I_{ays}, \tau_{ays}^2 ;\ \log \hat{I}_{ays} \Big)$$

giving the total log-likelihood

$$\ell = \ell_C + \ell_I$$

which is defined in terms of the strictly positive quantites, $M_{ay}$, $F_{ay}$, $Q_{ays}$ and $R_{y}$, and the observation variances $\sigma_{ay}$ and $\tau_{ays}$. As such, the log-likelihood is over-parameterised as there are many more parameters than observations. In order to reduce the number of parameters, $M_{ay}$ is assumed known (as is common). 

%====================================================================
%THE FOLLOWING NEEDS REVISION, need to bring in N1 submod
%====================================================================

The remaining parameters are written in terms of a linear combination of covariates $x_{ayk}$, e.g.

$$\log F_{ay} = \sum_k \beta_k x_{ayk}$$

where $k$ is the number of parameters to be estimated and is sufficiently small. Using this tecnique the quantities $\log F$, $\log Q$, $\log \sigma$ and $\log \tau$
%$\log \text{initial\,age\,structure}$ % this is not present in the above
(in bold in the equations above) can be described by a reduced number of parameters. The following section has more discussion on the use of linear models in \aFa.

%====================================================================

The \aFa statistical catch-at-age model can addionally allow for a functional relationship that links predicted recruitment $\tilde{R}$ based on spawning stock biomass and modelled recruitment $R$, to be imposed as a fixed variance random effect. [NEEDS REVISION, sentence not clear]

Options for the relationship are the hard coded models Ricker, Beverton Holt, smooth hockeystick or geometric mean. This is implemented by including a third component in the log-likelihood:

$$\ell_{SR} = \sum_y \ell_N \Big( \log \tilde{R}_y(a, b), \phi_y^2 ;\ \log R_y \Big)$$

giving the total log-likelihood

$$\ell = \ell_C + \ell_I + \ell_{SR}$$

Using the (time varying) Ricker model as an example, predicted recruitment is

$$\tilde{R}_y(a_y,b_y) = a_y S_{y-1} e^{-b_y S_{y-1}}$$

where $S$ is spawning stock biomass derived from the model parameters $F$ and $R$, and the fixed quantites $M$ and mean weights by year and age. It is assumed that $R$ is log-normally distributed, or equivalently, normally distributed on the log-scale about the (log) recruitment predicted by the SR model $\tilde{R}$, with known variance $\phi^2$, i.e.

$$\log R_y \sim \text{Normal} \Big( \log \tilde{R}_y, \phi_y^2 \Big)$$

which leads to the definition of $\ell_{SR}$ given above. In all cases $a$ and $b$ are strictly positive, and with the quantities $F$, $R$, etc. linear models are used to parameterise $\log a$ and/or $\log b$, where relevant.

%====================================================================
%THE FOLLOWING NEEDS REVISION, this is not just the default I guess, 
% it's always present since R predictions will be a mix of S/R and $\gamma$
%====================================================================

By default, recruitment $R$ as apposed to the reruitment predicted from a stock recruitment model $\tilde{R}$, is specified as a linear model with a parameter for each year, i.e.

$$\log R_y = \gamma_y$$

This is to allow modelled recruitment $R_y$ to be shrunk towards the stock recruitment model. However, if it is considered appropriate that recruitment can be determined exactly by a relationship with covariates, it is possible, to instead define $\log R$ in terms of a linear model in the same way as $\log F$, $\log Q$, $\log \sigma$ and $\log \tau$.  %But this is pretty much the same as taking a geometric mean, with a model on log a, and making the variance very small.

%====================================================================
% WE NEED TO ADD SOMETHING ABOUT HOW THE PLUSGROUP IS MODELLED
%====================================================================

%====================================================================

\section{Submodel structure}\label{sec:submod}

The \aFa stock assessment framework allows the user to set up a large number of different models. The mechanics which provide this flexibility are designed around the concept of submodels. Each unknown variable that must be estimated is treated as a linear model, for which the user has to define the model structure using \R formulas, including \href{http://cran.r-project.org/web/packages/mgcv/index.html}{\code{mgcv}} gam formulas. All submodels use the same specification process, the \R formula interface, wich gives lot's of flexibility to explore models and combination of submodels.

There are 5 submodels in operation:
\begin{itemize}
  \item a model for F-at-age ($F_{ay}$),
  \item a (list) of model(s) for abundance indices catchability-at-age ($Q_{ays}$),
  \item a model for recruitment ($R_y$),
  \item a list of models for the observation variance of catch-at-age and abundance indices ($\{\sigma^2_{ay}, \tau^2_{ays}\}$),
  \item a model for the initial age structure $N_{a,y=1}$,
\end{itemize}

When setting the structure of each submodel the user is in fact building the predictive model and its parameters. The optimization process, done through \ADMB, estimates the parameters and their variance-covariance matrix, allowing further analysis to be carried out, like simulation, prediction, diagnostics, etc. All the statistical machinery will be at the user's reach.

\subsection{Submodel building blocks and fundamental R formulas}

The elements available to build submodels formulas are 'age' and 'year', which can be used to build models with different structures. 

In \R's linear modelling language, a constant model is coded as \code{\midtilde 1}, while a slope over time would simply be \code{\midtilde year}, a smoother over time \code{\midtilde s(year, k=10)}, a model with a coefficient for each year (also called dummy variables) would be \code{\midtilde factor(year)}. Transformations of the variables are as usual, e.g. \code{\midtilde sqrt(year)}, etc; while combinations of all the above can be done non-convergence will limit the possibilities. 

Using the $F$ submodel as example the following specifies the models described in the previous paragraph:

<<fund_forms, fig.cap="Example of fundamental R formulas">>=
# models
m1 <- ~1
m2 <- ~ year
m3 <- ~ s(year, k=10)
m4 <- ~ factor(year)
m5 <- ~ sqrt(year)

# fits
fit1 <- sca(ple4, ple4.indices, fmodel=m1, fit="MP")
fit2 <- sca(ple4, ple4.indices, fmodel=m2, fit="MP")
fit3 <- sca(ple4, ple4.indices, fmodel=m3, fit="MP")
fit4 <- sca(ple4, ple4.indices, fmodel=m4, fit="MP")
fit5 <- sca(ple4, ple4.indices, fmodel=m5, fit="MP")

# plot
lst <- FLStocks(constant=ple4+fit1, linear=ple4+fit2, smooth=ple4+fit3, factor=ple4+fit4, sqrt=ple4+fit5)
lst <- lapply(lst, fbar)
lgnd <- list(points=FALSE, lines=TRUE, space='right')
xyplot(data~year, groups=qname, lst, auto.key=lgnd, type='l', ylab='fishing mortality')
@

The models above and their combinations can be used to model both 'age' and 'year'. The corresponding fits for age are:

<<fund_forms_age, fig.cap="Example of fundamental R formulas">>=
# models
m1 <- ~1
m2 <- ~ age
m3 <- ~ s(age, k=3)
m4 <- ~ factor(age)
m5 <- ~ sqrt(age)

# fits
fit1 <- sca(ple4, ple4.indices, fmodel=m1, fit="MP")
fit2 <- sca(ple4, ple4.indices, fmodel=m2, fit="MP")
fit3 <- sca(ple4, ple4.indices, fmodel=m3, fit="MP")
fit4 <- sca(ple4, ple4.indices, fmodel=m4, fit="MP")
fit5 <- sca(ple4, ple4.indices, fmodel=m5, fit="MP")

# plot
lst <- FLStocks(constant=ple4+fit1, linear=ple4+fit2, smooth=ple4+fit3, factor=ple4+fit4, sqrt=ple4+fit5)
lst <- lapply(lst, function(x) harvest(x)[,'2000'])
xyplot(data~age, groups=qname, lst, auto.key=lgnd, type='l', ylab='fishing mortality in 2000')
@

\subsection{The major effects available for modelling}

Although the building blocks for formulas are 'age' and 'year', in fact there are three effects that can be modelled for each submodel: 'age', 'year' and 'cohort'. As examples note the following models for fishing mortality.

<<>>=
# the age effect
ageeffect <- ~ factor(age)

# the year effect
yeareffect <- ~ factor(year)

# the cohort
cohorteffect <- ~ factor(year-age)

# the fits
fit1 <- sca(ple4, ple4.indices, fmodel=yeareffect)
fit2 <- sca(ple4, ple4.indices, fmodel=ageeffect)
fit3 <- sca(ple4, ple4.indices, fmodel=cohorteffect)
@

and the graphical representation of the three models in Figures~\ref{fig:majeffy} to \ref{fig:majeffc}.

<<majeffy, fig.cap="Major effects: the year effect (~ factor(year))">>=
wireframe(harvest(fit1), main='year effect')
@

<<majeffa, fig.cap="Major effects: the age effect (~ factor(age))">>=
wireframe(harvest(fit2), main='age effect')
@

<<majeffc, fig.cap="Major effects: the cohort effect (~ factor(year-age))">>=
wireframe(harvest(fit3), main='cohort effect')
@

\subsection{The submodel class and methods}

%====================================================================
% COLIN TO CHECK THE SECTION
%====================================================================
%
% Although the specification of each submodel is done through a R formula, internally the \aFa sca fit creates an object (of class 'submodel') which stores more information and allows a wide number of methods to be applied. The most important cases are prediction and simulation methods.
%
% The submodel class is a S4 class with the following slots:
%
% <<>>=
% showClass("submodel")
% @
%
% Objects of class 'submodel' are created in the fitting process using the formulas set by the user (or defaults if missing) and information from the fit. 
%
% For example begining with a simple \aFa fit:
%
% <<>>=
% # fit a model with indices "IBTS_Q1" and "IBTS_Q3"
% fit0 <- sca(ple4, ple4.indices[c("IBTS_Q1", "IBTS_Q3")],
%             fmodel = ~ s(age, k = 5),
%             qmodel = list( ~ s(age, k = 4), ~ s(age, k = 4)),
%             srmodel = ~ 1,
%             n1model = ~ s(age, k = 5),
%             vmodel = list( ~ 1, ~ 1, ~ 1),
%             verbose = FALSE)
% fit0
% @
%
% Within the sca fit object there's a collection of submodels for each component of the \aFa stock assessment model (fishing mortality, survey catchability, stock recruitment relationship, initial population and the observation variance for catch numbers and survey indices). For readability the following example will use the fishing mortality submodel although the methods and models described can be applied to each of the 5 submodels.
%
% The submodel can be extracted using:
%
% <<>>=
% fmod <- fmodel(fit0)
% @
%
% inside this object is a formula, the coeficients estimated by the model and their covariance matrix.
%
% <<sim_pred_fmodel_coef, results = 'hide', eval=FALSE>>=
% coef(fmod)
% vcov(fmod)
% @
%
% it is also possible to get the data behind the fit and the design matrix required to get the fitted values
%
% <<results = 'hide', eval=FALSE>>=
% as.data.frame(fmod)
% getX(fmod)
% @
%
% this makes it possble to check the fit manually like this (note that the data is centered in the model, which is addressed by adding the centering after the prediction is made, see code below):
%
% <<eval=FALSE>>= 
% dat <- as.data.frame(fmod)
% X <- getX(fmod)
% b <- coef(fmod)
% dat$fit <- exp(c(X %*% b))# + dat$data.centering)
%
% plot(fit ~ age, type = "l", data = dat[dat$year == 2016,],
%      ylab = "Estimated fishing mortality at age",
%      ylim = c(0, max(dat$fit)), las = 1)
% @
%
% or simulate from the fit by resampling based on the variance matrix of the coefficients, which is done by simulating many coefficients (here called `bsim`).
%
% <<eval=FALSE>>=
% # get the variance matrix of the coefficients and simulate
% Sigma <- vcov(fmod)[,,1]
% bsim <- mvrnorm(999, b, Sigma)
%
% fit_sim <- exp(X %*% bsim)
% dat$ci_lower <- apply(fit_sim, 1, quantile, p = 0.025)
% dat$ci_upper <- apply(fit_sim, 1, quantile, p = 0.975)
%
% plot(fit ~ age, type = "l", data = dat[dat$year == 2016,],
%      ylab = "Estimated fishing mortality at age",
%      ylim = c(0, max(dat$ci_upper)), las = 1)
% lines(ci_lower ~ age, data = dat[dat$year == 2016,], lty = 2)
% lines(ci_upper ~ age, data = dat[dat$year == 2016,], lty = 2)
% @
%
% The above code is to show how it is possible to use an FLa4a submodel to predict and simulate, and not intended to be used by the user.  The recomended way to predict from a submodel and make simulations is to use the `genFLQuant` function (and for finer control, there is also a `simulate` function).  The above example can be done using FLa4a functions as follows:
%
% <<eval=FALSE>>=
% fmod_fit <- genFLQuant(fmod)
% dat <- as.data.frame(fmod_fit)
%
% plot(data ~ age, type = "l", data = dat[dat$year == 2016,],
%      ylab = "Estimated fishing mortality at age",
%      ylim = c(0, max(dat$data)), las = 1)
% @
%
% or using the `ggplotFL` package:
%
% <<eval=FALSE>>=
% fmod_fit <- genFLQuant(fmod)
%
% ggplot(data = fmod_fit[,"2016"], aes(x = age, y = data)) + 
%   geom_point() + geom_line() + 
%   ylab("Estimated fishing mortality at age") +
%   ylim(0, max(fmod_fit))
% @
%
% and the simulations and confidence intervals are easily computed:
%
% <<eval=FALSE>>=
% # simulate 999 
% fmod_fit_sim <- genFLQuant(fmod, nsim = 999)
%
% # reduce to quantiles
% fmod_fit_sim <- quantile(fmod_fit_sim[,"2016"], c(0.025, 0.50, 0.975))
% @
%
% plotting can be done by converting to a data.frame, reshaping and using standard `ggplot2` functionality
%
% <<eval=FALSE>>=
% dat <- 
%   reshape(
%     as.data.frame(fmod_fit_sim, drop=TRUE), 
%     timevar = "iter", idvar = "age", direction = "wide"  
%   )
%
% ggplot(data=dat, aes(x = age, y = `data.50%`)) +
%   geom_ribbon(aes(ymin = `data.2.5%`, ymax = `data.97.5%`), 
%               fill="red", alpha = .15) +
%   geom_point() + geom_line() + 
%   ylab("Estimated fishing mortality at age") +
%   ylim(0, max(fmod_fit_sim))
% @
%
\section{The statistical catch-at-age stock assessment framework}\label{sec:sca}

The \aFa stock assessment framework is implemented in \R through the method \code{sca()}. The method call requires as a minimum a FLStock object and a FLIndices object, in which case the default submodels will be set by the method.

%Having described building blocks, basic formulations and effects available to build a submodel's model, it's important to look into specific formulations and relate them to commonly known representations. Note that although a large number of formulations are available for each submodel, the user must carefuly decide on the full stock assessment model being build and avoid over-paramerizing. Over-parametrization may lead to non-convergence, but may also end up not being very useful for prediction/forecasting, which is one of the main objectives of stock assessment.

<<>>=
data(ple4)
data(ple4.indices)
fit <- sca(ple4, ple4.indices)
stk <- ple4 + fit
plot(stk)
@

By calling the fitted object the default submodel formulas are printed in the console:

<<>>=
fit
@

To set specific submodels the user has to write the relevant R formula and include it in the call. The arguments for each submodel are self-explanatory: fishing mortality is 'fmodel', indices' catchability is 'qmodel', stock-recruitment is 'srmodel', observation variance is 'vmodel' and for initial year's abundance is 'n1model'. The following model comes closer to the official stock assessment of North Sea plaice, as such we'll name it $0$ and keep it for future comparisons:


For future referencing we'll start with a base fit to be used for future comparisons, named fit 0.

<<fit0>>=
@


<<>>=
fmod0 <- ~s(age, k=6)+s(year, k=10)+te(age, year, k=c(3,8))
qmod0 <- list(~s(age, k = 4), ~s(age, k = 3), ~s(age, k = 3)+year, ~s(age, k = 3), ~s(age, k = 4), ~s(age, k = 6))
srmod0 <- ~ s(year, k=20)
vmod0 <- list(~s(age, k=4), ~1,  ~1, ~1, ~1, ~1, ~1, ~1)
n1mod0 <- ~ s(age, k=3)
fit0 <- sca(ple4, ple4.indices, fmodel=fmod0, qmodel=qmod0, srmodel=srmod0, n1model=n1mod0, vmodel=vmod0)
stk0 <- ple4 + fit0
plot(stk0)
@

As before by calling the fitted object submodels' formulas are printed in the console:

<<>>=
fit
@

The method \code{sca} has other arguments which may be set by the user:

\begin{description}
	\item [covar:] a \code{FLQuant} with covariates; 
	\item [wkdir:] a folder (character) where the \ADMB files will be saved for posterior inspection by the user;
	\item [verbose:] be more verbose (logical);
	\item [fit:] type of fit (character), 
	\begin{itemize}
	  \item 'MP' runs the minimizer without trying to invert the hessian and as such doesn't return the covariance matrix of the parameters, normally used inside \MSE loops where parameter variance may not be relevant; 
	  \item 'assessment' runs minimizer and inverts hessian, returns the covariance matrix of the estimated parameters and the convergence criteria set in \ADMB; 
	  \item 'MCMC' runs \ADMB's MCMC fit
	\end{itemize} 
	\item [center:] shall observations be centered before fitting (logical);
	\item [mcmc:] \ADMB's MCMC arguments (character vector), must be paired with \code{fit="MCMC"}. 
\end{description}

There are a set of methods for \aFa fit objects which help manipulating \code{sca()} results, namely:

\begin{description}
	\item [+:] update the stock object with the fitted fishing mortalities, population abundance and catch in numbers at age; 
\end{description}


\subsection{Fishing mortality submodel ($F_{ay}$)}

<<echo=FALSE>>=
data(ple4)
data(ple4.indices)
@

We will now take a look at some examples for $F$ models and the forms that we can get.

% A non-separable model, where we consider age and year to interact can be modeled using a smooth interaction term in the F model using a tensor product of cubic splines with the `te` method (`r fign('te1')`), again borrowed from [mgcv](http://cran.r-project.org/web/packages/mgcv/index.html).
%
% <<>>=
% fmod <- ~ te(age, year, k = c(4,20))
% fit <- sca(ple4, ple4.indices[1], fmod)
% @
%
% <<te1, fig.cap="Fishing mortality smoothed non-separable model">>=
% wireframe(harvest(fit), zlab="F")
% @
%
% In the last examples the fishing mortalities (Fs') are linked across age and time.  What if we want to free up a specific age class because in the residuals we see a consistent pattern.  This can happen, for example, if the spatial distribution of juveniles is disconnected to the distribution of adults.  The fishery focuses on the adult fish, and therefore the the F on young fish is a function of the distribution of the juveniles and could deserve a specific model. This can be achieved by adding a component for the year effect on age 1 (`r fign('age1')`).
%
% <<>>=
% fmod <- ~ te(age, year, k = c(4,20)) + s(year, k = 5, by = as.numeric(age==1))
% fit <- sca(ple4, ple4.indices[1], fmod)
% @
%
% <<age1, fig.cap="Fishing mortality age-year interaction model with extra age 1 smoother.">>= 
% wireframe(harvest(fit), zlab="F")
% @
%
\subsubsection{Separable model}

One of the most useful models for fishing mortality is one in which 'age' and 'year' effects are independent, that is, where the shape of the selection pattern does not change over time, but the overall level of fishing mortality do. Commonly called a 'separable model'. 

A full separable model in \aFa is written using the \code{factor} function which converts age and year effects into categorical values, forcing a different coefficient to be estimated for each level of both effects. This model has \code{age x year} number of parameters.

<<>>=
fmod1 <- ~ factor(age) + factor(year)
fit1 <- sca(ple4, ple4.indices, fmodel=fmod1, fit="MP")
@

One can reduce the number of parameters and add dependency along both effects, although still keeping independence of each other, by using smoothers rather than \code{factor}. We'll use a (unpenalised) thin plate spline provided by package \href{http://cran.r-project.org/web/packages/mgcv/}{mgcv} method \code{s()}. We're using the North Sea Plaice data, and since it has 10 ages we will use a simple rule of thumb that the spline should have fewer than $\frac{10}{2} = 5$ degrees of freedom, and so we opt for 4 degrees of freedom. We will also do the same for year and model the change in $F$ through time as a smoother with 20 degrees of freedom.

<<>>=
fmod2 <- ~ s(age, k=4) + s(year, k=20)
fit2 <- sca(ple4, ple4.indices, fmodel=fmod2, fit="MP")
@

An interesting extension of the separable model is the 'double separable' where a third factor or smoother is added for the cohort effect.

<<>>=
fmod3 <- ~ s(age, k=4) + s(year, k=20) + s(as.numeric(year-age), k=10)
fit3 <- sca(ple4, ple4.indices, fmodel=fmod3, fit="MP")
@

Figures~\ref{fig:sep00} and \ref{fig:sep01} depicts the three models selectivities for each year. Each separable model has a single selectivity that changes it's overall scale in each year, while the double separable introduces some variability over time by modeling the cohort factor.

<<sep00, fig.cap="Selection pattern of separable models. Each line represents the selection pattern in a specific year. Independent age and year effects (factor), internally dependent age and year (smooth), double separable (double).">>=
flqs <- FLQuants(factor=harvest(fit1), smooth=harvest(fit2), double=harvest(fit3))
pset <- list(strip.background=list(col="gray90"))
xyplot(data~age|qname, groups=year, data=flqs, type="l", col=1, layout=c(3,1), ylab="fishing mortality", par.settings=pset)
@

<<sep01, fig.cap="Fishing mortality of separable models. Independent age and year effects (factor), internally dependent age and year (smooth), double separable (double).">>=
wireframe(data~age+year|qname, data=as.data.frame(flqs), layout=c(3,1))
@

\subsubsection{Constant selectivity for contiguous ages or years}

To set these models we'll use the method \code{replace()} to define which ages or years will be modelled together with a single coefficient. The following example shows \code{replace()} in operation. The dependent variables used in the model will be changed and attributed the same age or year, as such during the fit observations of those age or year with will be seen as replicates. One can think of it as sharing the same mean value, which will be estimated by the model.

<<>>=
age <- 1:10
# last age same as previous
replace(age, age>9, 9)
# all ages after age 6
replace(age, age>6, 6)

year <- 1950:2010
replace(year, year>2005, 2005)
@

In the $F$ submodel one can use this method to fix the estimation of $F$ in the plus group to be the same as in the last non-aggregated age.

<<>>=
fmod <- ~ s(replace(age, age>9, 9), k=4) + s(year, k=20)
fit <- sca(ple4, ple4.indices, fmod)
@

<<ctsselage, fig.cap="F-at-age fixed above age 9">>=
wireframe(harvest(fit), zlab="F")
@

Or estimate the average $F$ in the most recent years, instead of averaging after the assessment to compute the \emph{statu quo} selection pattern.

<<>>=
fmod <- ~ s(age, k=4) + s(replace(year, year>2013, 2013), k=20)
fit <- sca(ple4, ple4.indices, fmod)
@

<<ctsselyear, fig.cap="F-at-age fixed for the most recent 5 years">>=
wireframe(data~age+year, data=harvest(fit)[,ac(2010:2017)], screen=c(z=-130, y=0, x=-60), zlab="F")
@

\subsubsection{Time blocks selectivity}

To define blocks of data \code(sca()) uses the method \code{breakpts()}, which creates a factor from a vector with levels defined by the second argument.

<<>>=
year <- 1950:2010
# two levels separated in 2000
breakpts(year, 2000)
# five periods with equal interval
breakpts(year, seq(1949, 2010, length=6))
@ 

Note \code{seq()} computes 'left-open' intervals, which means that to include 1950 the sequence must start one year earlier. 

These methods can be used to create discrete time series, for which a different selection pattern is allowed in each block. This is called an interaction in statistical modelling parlance, and typically a \code{*} denotes an interaction term; for smoothers an interaction is achieved using the \code{by} argument. When this argument is a \code{factor} a replicate of the smooth is produced for each factor level. 

In the next case we'll use the \code{breakpts()} to split the time series at 1990, although keeping the same shape in both periods, a thin plate spline with 3 knots (Figure~\ref{fig:brk}).

<<>>=
fmod <- ~s(age, k = 3, by = breakpts(year, 1990))
fit <- sca(ple4, ple4.indices, fmod)
@

<<brk, echo=FALSE, fig.cap="F-at-age in two periods using in both cases a thin plate spline with 3 knots">>=
wireframe(harvest(fit), zlab="F")
@

\subsubsection{Time changing selectivity}

In many cases, it may be desirable to allow the selection pattern to evolve over time, from year to year. Again there are several ways to do this, one way is to estimate a mean selection pattern, while also allowing F to vary over time for each age. This is like a seperate smoother over year, with 'age blocks' so, looking back at previous examples, we have:

<<>>= 
fmodel <- ~ s(year, k = 15, by = factor(age)) + s(age, k = 4)
@

This is a type of interaction between age and year, but the only connection (or correlation) across ages is via the smoother on age, however there are still 15 degrees of freedom for each age, so the model 5 x 15 + 4 = 69 degrees of freedom.  To include correlation across ages and years together then the tensor product (`te()` function) is used, this has the effect of restricting the flexibility of the model for F.  In the following, there is a smoother in 2 dimensions (age and year) where there is 5 degrees of freedom in the age direction, and 15 in the year dimension, resulting in a total of 5 x 15 = 65 degrees of freedom

<<>>=
fmodel <- ~ te(age, year, k = c(5, 15))
@

Often the above formulations provide too much flexibility, and a more complicated, but simpler model is preferable:

<<>>=
fmodel <- ~ s(age, k = 4) + s(year, k = 15) + te(age, year, k = c(3, 5))
@

in the above model, the main effects for age and year still have similar flexibility to the full tensor model, however, the interaction (or the change in F at age over time) has been restricted, so that the full model now has 4 + 15 + 3 x 5 = 34 degrees of freedom.

\subsubsection{Trawl fleets}

\subsubsection{Nets and Liners fleets}

\subsubsection{Multigear fleets}

\subsubsection{Trawl surveys}

\subsubsection{Closed form selection pattern}

One can use a closed form for the selection pattern. The only requirement is to be able to write it as a \R formula, the example below uses a logistic form.

<<>>=
fmod <- ~ I(1/(1+exp(-age)))
fit <- sca(ple4, ple4.indices, fmod)
@

<<logistic, fig.cap="F-at-age logistic">>=
wireframe(harvest(fit), zlab="F")
@

% \subsubsection{More models}
%
% More complicated models can be built with these tools. For example, `r fign('ageind')` shows a model where the age effect is modelled as a smoother (the same thin plate spline) throughout years but independent from each other.
%
% <<>>=
% fmod <- ~ factor(age) + s(year, k=10, by = breakpts(age, c(2:8)))
% fit <- sca(ple4, ple4.indices, fmod)
% @
%
% <<ageind, fig.cap="F-at-age as thin plate spline with 3 knots for each age">>=
% wireframe(harvest(fit), zlab="F")
% @
%
% A quite complex model that implements a cohort effect can be set through the following formula. `r fign('coh')` shows the resulting fishing mortality. Note that in this case we end up with a variable F pattern over time, but rather than using 4 * 10 = 40 parameters, it uses, 4 + 10 + 10 = 24.
%
% <<>>=
% fmodel <- ~ s(age, k = 4) + s(pmax(year - age, 1957), k = 10) + s(year, k = 10)
% fit <- sca(ple4, ple4.indices, fmodel=fmodel)
% @
%
% <<coh, echo=FALSE, fig.cap="F-at-age with a cohort effect.">>=
% wireframe(harvest(fit), zlab="F")
% @
%
\subsection{Abundance indices catchability submodel ($Q_{ays}$)}

The catchability submodel is set up the same way as the $F$ submodel and the tools available are the same. The only difference is that the submodel is set up as a list of formulas, where each formula relates with one abundance index. There's no limitation in the number of indices or type that can be used for a fit. It's the analyst that has to decide based on her/his expertise and knowledge of the stock and fleet dynamics.

\subsubsection{Catchability submodel for age based indices}

A first model is simply a dummy effect on age, which means that a coefficient will be estimated for each age. Note that this kind of model considers that levels of the factor are independent (Figure~\ref{fig:dummyage}).

<<dummyage, fig.cap="Catchability age independent model">>=
qmod <- list(~factor(age))
fit <- sca(ple4, ple4.indices[1], qmodel=qmod)
qhat <- predict(fit)$qmodel[[1]]
wireframe(qhat, zlab="q")
@

If one considers catchability at a specific age to be dependent on catchability on the other ages, similar to a selectivity modelling approach, one option is to use a smoother at age, and let the data 'speak' regarding the shape (Figure~\ref{fig:smoothage}).

<<smoothage, fig.cap="Catchability smoother age model">>=
qmod <- list(~ s(age, k=4))
fit <- sca(ple4, ple4.indices[1], qmodel=qmod)
qhat <- predict(fit)$qmodel[[1]]
wireframe(qhat, zlab="q")
@

Finally, one may want to investigate a trend in catchability with time, very common in indices built from CPUE data. In the example given here we'll use a linear trend in time, set up by a simple linear model (Figure~\ref{fig:qtrend}).

<<qtrend, fig.cap="Catchability with a linear trend in year">>=
qmod <- list( ~ s(age, k=4) + year)
fit <- sca(ple4, ple4.indices[1], qmodel=qmod)
qhat <- predict(fit)$qmodel[[1]]
wireframe(qhat, zlab="q")
@

\subsubsection{Catchability submodel for age aggregated biomass indices}

The previous section was focused on age disaggregated indices, but age aggregated indices (CPUE, biomass, DEPM, etc) may also be used to tune the total biomass of the population. In these cases a different class for the index must be used, the \code{FLIndexBiomass}, which uses a vector \code{index} with the age dimension called 'all'. Note that in this case the qmodel should be set without age factors, although it can have a 'year' component and covariates if needed. An interesting feature with biomass indices is the age range they refer to can be specified.

<<>>=
# simulating a biomass index (note the name of the first dimension element) using 
# the ple4 biomass and an arbritary catchability of 0.001 plus a lognormal error.
dnms <- list(age="all", year=range(ple4)["minyear"]:range(ple4)["maxyear"])
bioidx <- FLIndexBiomass(FLQuant(NA, dimnames=dnms))
index(bioidx) <- stock(ple4)*0.001
index(bioidx) <- index(bioidx)*exp(rnorm(index(bioidx), sd=0.1))
range(bioidx)[c("startf","endf")] <- c(0,0)
# note the name of the first dimension element
index(bioidx)
# fitting the model
fit <- sca(ple4, FLIndices(bioidx), qmodel=list(~1))
@

To estimate a constant selectivity over time one used the model \code{\midtilde 1}. As a matter of fact the estimate value, \Sexpr{round(predict(fit)$qmodel[[1]][1,drop=TRUE],5)}, is not very far from the simulated one, 0.001.

An example where the biomass index refers only to age 2 to 4 (for example a CPUE that targets these particular ages).

<<>>=
# creating the index
dnms <- list(age="all", year=range(ple4)["minyear"]:range(ple4)["maxyear"])
bioidx <- FLIndexBiomass(FLQuant(NA, dimnames=dnms))
# but now use only ages 2:4
index(bioidx) <- tsb(ple4[ac(2:4)])*0.001
index(bioidx) <- index(bioidx)*exp(rnorm(index(bioidx), sd=0.1))
range(bioidx)[c("startf","endf")] <- c(0,0)
# to pass this information to the model one needs to specify an age range
range(bioidx)[c("min","max")] <- c(2,4)
# fitting the model
fit <- sca(ple4, FLIndices(bioidx), qmodel=list(~1))
@

Once more the estimate value, \Sexpr{round(predict(fit)$qmodel[[1]][1,drop=TRUE],5)}, is not very far from the simulated one, 0.001.

\subsubsection{Catchability submodel for single age indices}

Similar to age aggregated indices one may have an index that relates only to one age, like a recruitment index. In this case the \code{FLIndex} object must have in the first dimension the age it referes to. The fit is then done relating the index with the proper age in numbers. Note that in this case the qmodel should be set without age factors, although it can have a 'year' component and covariates if needed.

<<>>=
idx <- ple4.indices[[1]][1]
fit <- sca(ple4, FLIndices(recidx=idx), qmodel=list(~1))
# the estimated catchability is
predict(fit)$qmodel[[1]]
@

\subsection{Stock-recruitment submodel ($R_y$)}

The S/R submodel is a special case, in the sense that it can be set up with the same linear tools as the $F$ and $Q$ models, but it can also use some hard coded models. The example shows how to set up a simple dummy model with \code{factor()}, a smooth model with \code{s()}, a Ricker model (\code{ricker()}), a Beverton and Holt model (\code{bevholt()}), a hockey stick model (\code{hockey()}), and a geometric mean model (\code{geomean()}). See Figure~\ref{fig:srmod} for results. As mentioned before, the 'structural' models have a fixed variance, which must be set by defining the coefficient of variation.

<<>>=
srmod <- ~ factor(year)
fit <- sca(ple4, ple4.indices, srmodel=srmod)
srmod <- ~ s(year, k=10)
fit1 <- sca(ple4, ple4.indices, srmodel=srmod)
srmod <- ~ ricker(CV=0.05)
fit2 <- sca(ple4, ple4.indices, srmodel=srmod)
srmod <- ~ bevholt(CV=0.05)
fit3 <- sca(ple4, ple4.indices, srmodel=srmod)
srmod <- ~ hockey(CV=0.05)
fit4 <- sca(ple4, ple4.indices, srmodel=srmod)
srmod <- ~ geomean(CV=0.05)
fit5 <- sca(ple4, ple4.indices, srmodel=srmod)
@

<<srmod, fig.cap="Stock-recruitment models fits">>=
flqs <- FLQuants(factor=stock.n(fit)[1], smother=stock.n(fit1)[1], ricker=stock.n(fit2)[1], bevholt=stock.n(fit3)[1], hockey=stock.n(fit4)[1], geomean=stock.n(fit5)[1])
xyplot(data~year, groups=qname, data=flqs, type="l", auto.key=list(points=FALSE, lines=TRUE, columns=3), ylab="No. recruits")
@

\subsection{Observation variance submodel ($\{\sigma^2_{ay}, \tau^2_{ays}\}$)}

The variance model allows the user to set up the shape of the observation variances $\sigma^2_{ay}$ and $\tau^2_{ays}$. This is an important subject related with fisheries data used for input to stock assessment models. 

The defaults assume a U-shape model for catch-at-age and constant variance for abundance indices. The first relies on the fact that it's common to have more precision on the most represented ages and less precision on the less frequent ages which tend to  be the younger and older individuals. These sizes are less caught by the fleets and as such do not appear as often at the auction markets samples. With regards to the abundance indices, one assumes a scientific survey to have a well designed sampling scheme and protocols which keep observation error at similar levels across ages.

<<>>=
vmod <- list(~s(age, k=3), ~1)
fit1 <- sca(ple4, ple4.indices[1], vmodel=vmod)
vmod <- list(~s(age, k=3), ~s(age, k=3))
fit2 <- sca(ple4, ple4.indices[1], vmodel=vmod)
@

Variance estimated for the constant model is \Sexpr{round(predict(fit)$vmodel[[2]][1,drop=TRUE],3)} while for the U-shape model, fitted with a smoother, changes with ages (Figure~\ref{fig:vmod}).

<<vmod, fig.cap="Abundance index observation variance estimate">>=
wireframe(predict(fit2)$vmodel[[2]], zlab="variance")
@

Observation variance options have an impact in the final estimates of population abundance, which can be seen in Figure~\ref{fig:vmodimpact}.  

<<vmodimpact, fig.cap="Population estimates using two different variance models">>=
flqs <- FLQuants(smother=stock.n(fit1), factor=stock.n(fit2))
xyplot(data~year|age, groups=qname, data=flqs, type="l",
       scales=list(y=list(relation="free", draw=FALSE)),
       auto.key=list(points=FALSE, lines=TRUE, columns=2),
       par.settings=list(superpose.line=list(col=c("gray35", "black")),
       strip.background=list(col="gray90")), ylab="")
@

\subsection{Initial year abundance submodel ($N_{a,y=1}$)}

The submodel for the stock number at age in the first year of the time series is set up with the usual modelling tools (Figure~\ref{fig:ny1}). Beare in mind that the year effect does not make sense here since it refers to a single year, the first in the time series of data available. This model has its influence limited to the initial lower triangle of the population matrix, which in assessments with long time series doesn't make much difference. Nevertheless, when modelling stocks with short time series in relation to the number of ages present, it becomes more important and should be given proper attention.

<<>>=
n1mod <- ~s(age, k=3)
fit1 <- sca(ple4, ple4.indices, fmodel=fmod0, qmodel=qmod0, srmodel=srmod0, vmodel=vmod0, n1model=n1mod)
n1mod <- ~factor(age)
fit2 <- sca(ple4, ple4.indices, fmodel=fmod0, qmodel=qmod0, srmodel=srmod0, vmodel=vmod0, n1model=n1mod)
flqs <- FLQuants(smother=stock.n(fit1)[,1], factor=stock.n(fit2)[,1])
@

<<ny1, fig.cap="Nay=1 models">>=
pset <- list(superpose.line=list(col=c("gray50", "black"), lty=c(1,2)))
xyplot(data~age, groups=qname, data=flqs, type="l", auto.key=lgnd, par.settings=pset, ylab="")
@

The impact in the overall perspective of the stock status is depicted in Figure~\ref{fig:n1modimpact}. As time goes by the effect of this model vanishes and the fits become similar.  

<<n1modimpact, fig.cap="Population estimates using two different variance models">>=
flqs <- FLQuants(smother=stock.n(fit1), factor=stock.n(fit2))
pset$strip.background <- list(col="gray90")
scl <- list(y=list(relation="free", draw=FALSE))
xyplot(data~year|factor(age), groups=qname, data=flqs, type="l", scales=scl, auto.key=lgnd, par.settings=pset, ylab="")
@

\subsection{Data weigthing}

%====================================================================
% COLIN TO CHECK THE SECTION
%====================================================================

By default the likelihood components are not weighted and the contribution of each to the maximum likelihood depends on their own likelihood score. However, the user may change these weights by penalizing data points, the $w_{ays}$ in section~\ref{sec:maths}. The likelihood score of each data point will be multiplied by the normalized weights ($\sum w_{ays} = 1$). This is done by adding a variance matrix to the \code{catch.n} and \code{index.n} slots of the stock and index objects. The values should be given as coefficients of variation on the log scale, so that variance is $\log{({CV}^2 + 1)}$. Figures Figure~\ref{fig:likwgt} and \ref{fig:likwgtimpact} show the results of the two fits in the population abundance and stock summary.  

<<>>=
stk <- ple4
idx <- ple4.indices[1]
# cv of observed catches
varslt <- catch.n(stk)
varslt[] <- 0.4
catch.n(stk) <- FLQuantDistr(catch.n(stk), varslt)
# cv of observed indices
varslt <- index(idx[[1]])
varslt[] <- 0.1
index.var(idx[[1]]) <- varslt
# run
fit1 <- sca(stk, idx, fmodel=fmod0, qmodel=qmod0, srmodel=srmod0, vmodel=vmod0, n1model=n1mod0)
flqs <- FLQuants(nowgt=stock.n(fit0), extwgt=stock.n(fit1))


@

<<likwgt, fig.cap="Stock summary of distinct likelihood weightings">>=
xyplot(data~year|factor(age), groups=qname, data=flqs, type="l", scales=scl, auto.key=lgnd, par.settings=pset, ylab="")
@

<<likwgtimpact, fig.cap="Population estimates using two different variance models">>=
flsts <- FLStocks(nowgt=ple4+fit0, wgt=ple4 + fit1)
plot(flsts)
@

Note that by using a smaller CV for the index, one is increasing the contribution of the survey and penalizing catch at age, in relative terms. The ratio between likelihood scores of both fits show this effect with catch at age increasing by 2.3 while the index increases almost 8 fold.

<<>>=
fit0 <- sca(ple4, ple4.indices[1], fmodel=fmod0, qmodel=qmod0, srmodel=srmod0, vmodel=vmod0, n1model=n1mod0)
(fitSumm(fit1)/fitSumm(fit0))[c(2,8,9),]
@

\subsection{Working with covariates}

In linear model one can use covariates to explain part of the variance observed on the data that the 'core' model does not explain. The same can be done in the \aFa framework. The example below uses the North Atlantic Oscillation (NAO) index to model recruitment.

<<>>=
nao <- read.table("https://www.cpc.ncep.noaa.gov/products/precip/CWlink/pna/norm.nao.monthly.b5001.current.ascii.table", skip=1, fill=TRUE, na.strings="-99.90")
dnms <- list(quant="nao", year=1950:2024, unit="unique", season=1:12, area="unique")
nao <- FLQuant(unlist(nao[,-1]), dimnames=dnms, units="nao")
nao <- seasonMeans(trim(nao, year=dimnames(stock.n(ple4))$year))
@

First by simply assuming that the NAO index drives recruitment (`r fign('naor')`).

<<>>=
srmod <- ~ nao
fit2 <- sca(ple4, ple4.indices[1], qmodel=list(~s(age, k=4)), srmodel=srmod, covar=FLQuants(nao=nao))
flqs <- FLQuants(simple=stock.n(fit)[1], covar=stock.n(fit2)[1])
@
 
<<naor, echo=FALSE, fig.cap="Recruitment model with covariates. Using the NAO index as a recruitment index.">>=
xyplot(data~year, groups=qname, data=flqs, type="l",
       auto.key=list(points=FALSE, lines=TRUE, columns=2),
       par.settings=list(superpose.line=list(col=c("gray35", "black"), lty=c(2,1), lwd=c(1,1.5)),
       strip.background=list(col="gray90")), ylab="")
@

In a second model we're using the NAO index not to model recruitment directly but to model one of the parameters of the S/R function (`r fign('naor2')`).

<<>>=
srmod <- ~ ricker(a=~nao, CV=0.25)
fit3 <- sca(ple4, ple4.indices[1], qmodel=list(~s(age, k=4)), srmodel=srmod, covar=FLQuants(nao=nao))
flqs <- FLQuants(simple=stock.n(fit)[1], covar=stock.n(fit3)[1])
@

<<naor2, echo=FALSE, fig.cap="Recruitment model with covariates. Using the NAO index as a covariate for the stock-recruitment model parameters.">>=
xyplot(data~year, groups=qname, data=flqs, type="l",
       auto.key=list(points=FALSE, lines=TRUE, columns=2),
       par.settings=list(superpose.line=list(col=c("gray35", "black"), lty=c(2,1), lwd=c(1,1.5)),
       strip.background=list(col="gray90")), ylab="")
@

Note that covariates can be added to any submodel using the linear model capabilities of R.

\subsection{Assessing \ADMB files}

The framework gives access to the files produced to run the \ADMB fitting routine through the argument \code{wkdir}. When set up all the \ADMB files will be left in the directory. Note that the \ADMB tpl file is distributed with the \code{FLa4a}. One can get it from your \R library, under the folder \code{myRlib/FLa4a/admb/}.

<<eval=FALSE>>=
fit1 <- sca(ple4, ple4.indices, wkdir="fit1run")
@

\subsection{Missing observations in the catch matrix or index}

%====================================================================
% COLIN TO CHECK THE SECTION
%====================================================================

How are the data interpolated etc ...

\section{Diagnostics}\label{sec:diagn}

There's a large number of diagnostics that can be computed for a stock assessment model, the \aFa framework implements several analysis of residuals, visualizations and statistics that can be used to evaluate the fit quality and chose across multiple fits.

\subsection{Residuals}

Residuals are a ubiquos metrics to check quality of a fit. For \code{sca()} fits there are out-of-the-box methods to compute in the log scale, raw residuals (aka deviances), standardized residuals and pearson residuals. A set of plots to inspect residuals and evaluate fit quality and assumptions are implemented.

Consider $x_{ay}$ to be either a catch-at-age matrix ($C_{ay}$) or one abundance index ($I_{ay}$\footnote{For simplicity of notation we'll avoid the subscript $s$ in $I$, since we're referring to individual indices}) and $d$ to represent residuals.

Raw residuals are compute by $d_{ay} = \log{x_{ay}} - \log{\tilde{x}_{ay}}$ and have distribution $N(0,\upsilon^2_{a})$. Standardized residuals will be compute with $d^s_{ay} = \frac{d_{ay}}{\hat{\upsilon}^2_{a}}$ where $\hat{\upsilon}^2_{a} = (n-1)^{-1} \sum_y(d_{ay})^2$. Pearson residuals scale raw residuals by the estimates of $\sigma^2$ or $\tau^2$, as such $d^p_{ay} = \frac{d_{ay}}{\tilde{\upsilon}^2_{a}}$ where $\tilde{\upsilon}^2_{a} = \tilde{\sigma}^2_{a}$ for catches, or $\tilde{\upsilon}^2_{a} = \tilde{\tau}^2_{a}$ for each index of abundance.  

The \code{residuals()} method will compute these residuals and generate a object which can be plotted using a set of packed methods. The argument \code{type} will allow the user to chose which residuals will be computed. By default the method computes standardized residuals.

<<>>=
fit <- sca(ple4, ple4.indices)
d_s <- residuals(fit, ple4, ple4.indices)
@

Figure~\ref{fig:res} shows a scatterplot of standardized residuals with a smoother to guide (or mis-guide ...) your visual analysis. Note that the standardization should produce residuals with variance ~1, which means that most residual values should be between $\sim -2$ and $\sim 2$.

<<res, fig.cap="Standardized residuals for abundance indices and catch numbers (catch.n). Each panel is coded by age class, dots represent standardized residuals and lines a simple smoother.">>=
plot(d_s)
@

When plotting residuals by default the auxiliar line is a smoother. However it's possible to use other type of lines by setting the argument "auxline" in plot. The argument can take the values used by xyplot, which are (from panel.xyplot help page) one or more of the following: "p", "l", "h", "b", "o", "s", "S", "g", "r", "a", "smooth", and "spline". If type has more than one element, an attempt is made to combine the effect of each of the components. The behaviour if any of the first five are included is similar to the effect of the corresponding type in plot: "p" and "l" stand for points and lines respectively; "b" and "o" (for ‘overlay’) plot both; "h" draws vertical (or horizontal if horizontal = TRUE) line segments from the points to the origin. Types "s" and "S" are like "l" in the sense that they join consecutive points, but instead of being joined by a straight line, points are connected by a vertical and a horizontal segment forming a ‘step’, with the vertical segment coming first for "s", and the horizontal segment coming first for "S".  "g" adds a reference grid. Type "r" adds a linear regression line, "smooth" adds a loess fit, "spline" adds a cubic smoothing spline fit, and "a" draws line segments joining the average y value for each distinct x value. Figure~\ref{fig:resaux} shows a regression line over the residuals instead of the loess smooother.  

<<resaux, fig.cap="Standardized residuals for abundance indices and catch numbers (catch.n). Each panel is coded by age class, dots represent standardized residuals and lines a simple smoother.">>=
plot(d_s, auxline="r")
@


The common bubble plot (\code{bubble()}) are shown in Figure~\ref{fig:bub}. It shows the same information as Figure~\ref{fig:res} but in a multivariate perspective.

<<bub, fig.cap="Bubbles plot of standardized residuals for abundance indices and for catch numbers (catch.n).">>=
bubbles(d_s)
@

Figure~\ref{fig:qq} shows a quantile-quantile plot to assess how well standardized residuals match a normal distribution.

<<qq, fig.cap="Quantile-quantile plot of standardized residuals for abundance indices and catch numbers (catch.n). Each panel is coded by age class, dots represent standardized residuals and lines the normal distribution quantiles.">>=
qqmath(d_s)
@

Pearson residuals can be computed and plotted the same way as standardized residuals by setting \code{fit='pearson'} (Figure~\ref{fig:resp}). These residuals are particularly useful to evaluate the \code{vmodel}, departures from $0$ should be fixed by tweaking the \code{vmodel} equations. 

<<resp, fig.cap="Pearson residuals for abundance indices and catch numbers (catch.n). Each panel is coded by age class, dots represent standardized residuals and lines a simple smoother.">>=
d_p <- residuals(fit, ple4, ple4.indices, type='pearson')
plot(d_p)
@

Finally, the raw residuals are computed by setting \code{fit='deviances'} and plotted the same way as before (Figure~\ref{fig:resr}). These residuals are usefull to identify which data points are not well modelled, showing a large dispersion of the residuals and requiring more attention from the analyst.

<<resr, fig.cap="Raw residuals for abundance indices and catch numbers (catch.n). Each panel is coded by age class, dots represent standardized residuals and lines a simple smoother.">>=
d_r <- residuals(fit, ple4, ple4.indices, type='deviances')
plot(d_r)
@

\subsection{Predictive skill}

An important feature of stock assessment model fits is the capacity to predict, since one of the most important analysis done with these fits is forecasting future fishing opportunities under pre-defined conditions. The \aFa framework implements a visualization of the fit's predictive skill for both catch-at-age and abundance indices. These are generated by the method \code{plot()} with the fit object and a \code{FLStock} (Figure~\ref{fig:selplt}) or \code{FLIndices} (Figure~\ref{fig:idxplt}) object as arguments.

<<selplt, fig.cap="Predict and observed catch-at-age">>=
plot(fit, ple4)
@

<<idxplt, fig.cap="Predict and observed abundance-at-age">>=
plot(fit, ple4.indices)
@

\subsection{Aggreagted catch in weight}

Although a statistical catch-at-age model assumes errors in catch-at-age and, as such, errors in the total catch in weight, there's still interest to evaluate how close the model estimates are of the observed catch in weight\footnote{Some analysts believe this is the most important diagnostic since total catch should be trusted. Needless to say we don't agree and consider reported catch in weight one of the less reliable pieces of information available for stock assessment.}. The implementation of this diagnopstics is done through the method \code{computeCatchDiagnostics()}, which can be visualized with \code{plot()} (Figure~\ref{c_d}).   

<<catchdiag, fig.cap="Diagnostics for age aggregated catch in weight">>=
c_d <- computeCatchDiagnostics(fit, ple4)
plot(c_d)
@

\subsection{Fit summary, information and cross-validation metrics}

To get information about the likelihood fit the method \code{fitSumm()} can be used to report  number of parameters (\code{npar}), megative log-likelkihood (\code{nlogl}), \ADMB maximum gradient par (\code{maxgrad}), number of observations (\code{nobs}), generalized cross validation score (\code{gcv}), convergence flag (\code{convergence}) and acceptance rate (\code{accrate}) relevant for MCMC fits only. The second part refers to the likelihood value for each component.

<<>>=
fitSumm(fit)
@

Information criteria based metrics are reported with the methods: 

<<>>=
AIC(fit)
BIC(fit)
@

\subsection{The package a4adiags}

The package \code{a4adiags} contains some additional diagnostics based on the \textcolor{red}{reference}. Runs test checks weather the residuals are randomly distributed. A "run" is a sequence of the same sign residuals. Few runs indicate a trend or a correlation in the residuals while too many runs may suggest overfitting. 

The primary output of a runstest is a p-value where: a high p value $(p\leq 0.05)$ suggests that the residuals are randomly distributed, a low p value indicates a non-random pattern in the residuals.

<<>>=
library(a4adiags)
theme_set(theme_bw())
fit <- sca(mut09, mut09.idx, fmod = ~factor(age) + s(year, k = 8))
res <- residuals(fit, mut09, mut09.idx)
@

<<idxrunstest, fig.cap="Runstest for the abundance index">>=
plotRunstest(fit, mut09.idx, combine = F) + theme_bw() + facet_wrap(~age)
@

<<catchrunstest, fig.cap="Runstest for the catch by age">>=
plotRunstest(catch.n(mut09), catch.n(mut09 + fit), combine = F) + theme_bw() + facet_wrap(~age)
@

Green shading indicates no evidence $(p <  0.05)$ and red shading evidence $(p  >0.05)$ to reject the hypothesis of a randomly distributed time-series of residuals, respectively. The shaded (green/red) area spans three residual standard deviations to either side from zero, and the red points outside of the shading violate the ‘ $3\sigma$ limit’ for that series.

%Hindcast is used to assess the prediction skill of the model by removing a couple of years from the end of the time series and forecasting for that years, to assess how close the projected values are to the ones assessed by the model. The hindcast can estimate forecast bias by comparing the forecasted values to the reference model estimates. The `a4adiags` package use a different approach, it is based on a technique proposed by Kell \textcolor{red}{ref} named hindcast cross-validation where the forecasted values are compared with the observed ones. Additionally mean absolute squared error is computed, a statistic for evaluating the prediction skill. MASE basically compares the model prediction skill against a random walk, i.e against the predicted value of a random process based only on the previous year's observation.

%<<xcval, fig.cap="Hindcasting and MASE statistic. A MASE score > 1 indicates that the average model forecasts are worse than a random walk. Conversely, a MASE score of 0.5 indicates that the model forecasts twice as accurately as a naïve baseline prediction; thus, the model has prediction skill">>=
%xval <- a4ahcxval(mut09, FLIndices(mut09.idx), nyears = 5, nsq = 3, fmodel = ~factor(age) + s(year, k = 8))
%plotXval2(xval$indices) + ggtitle(paste0("Hindcast"))
%@

\section{Predict and simulate}\label{sec:predsim}

% To predict and simulate ``R`` uses the methods `predict()` and `simulate()`, which were implemented in ``FLa4a`` in the same fashion.
%
% <<pred0, eval=FALSE>>=
% fit <- sca(ple4, ple4.indices[1], fit="assessment")
% @
%
% \subsection{Basic functions}
%
% Simulation and prediction in FLa4a is based on three functions: `simulate` and `genFLQuant`. 
%
% \subsubsection{simulate()}
%
% Unlike the stats function `stats::simulate`, `FLa4a::simulate` will return the same object as it was passed.  For example, if you simulate from a FLa4a fit, you will get an FLa4a fit object back in which, the coefficients of the object are simulations from the model. Likewise, if you call simulate on a submodel object you will get back a submodel object inwhich the coefficients are simulations from the model.  Simulations are always done by generating random draws from a multivariate normal distribution with mean given by the coeffients of the model, and variance matrix given by the estimated covariance matrix of the coefficients (in practice this is a submatrix of the inverse of the hessian matrix).
%
% Simulate works on several classes from full fits right down to the individual model compents, so if `my\_fit` is a fitted \aFa model, then `simulate(my\_fit, nsim = 100)` will return a new fitted \aFa model where the model coefficients now have 100 iters and are drawn from the full variance matrix of the fitted model.  Similarly, `simulate(fmodel(my\_fit), nsim = 100)` will return a submodel with the same formula as the fishing mortality model as in the \aFa fit but where the coefficients are simulated from the variance matrix of the relavent parameters.  
%
% \subsubsection{genFLQuant()}
%
% This is a special function who's purpose is to return an `FLQuant` or `FLQuants`.  It essentially provides predictions from a model and provides in them as FLQuants of the correct dimensions.  `genFLQuant` also has an argument `nsim` which if set to a value greater than zero produces simulated predictions from the model based on simulations of the model coefficients.
%
% \subsection{submodels}
%
% In an `sca` fit individual submodel objects are often combined into a collection of `submodels`, for example the models for survey catchability are a collection of submodels
%
% <<sim_pred_qsubmodels>>=
% qmod <- qmodel(fit0)
% @
%
% which contains a submodel for each survey catchability. Now (almost) the same code can be run as before to plot the estimates with confidence intervals, the difference with submodels is all the results will be `FLQuants` so, `lapply` must be used to do computations on the predictions for each submodel seperately. 
%
% <<sim_pred_qmodel_predict_a4a, eval = FALSE>>=
% qmod_fit_sim <- genFLQuant(qmod, nsim = 999)
%
% # reduce to quantiles
% qmod_fit_sim <- lapply(qmod_fit_sim, "[", j = "2016")
% qmod_fit_sim <- lapply(qmod_fit_sim, quantile, prob = c(0.025, 0.50, 0.975))
%
% # reshape
% dat <- 
%   reshape(
%     as.data.frame(qmod_fit_sim, drop=TRUE), 
%     timevar = "iter", idvar = c("age", "qname"), direction = "wide"  
%   )
%
% # plot
% ggplot(data=dat, aes(x = age, y = `data.50%`)) +
%   geom_ribbon(aes(ymin = `data.2.5%`, ymax = `data.97.5%`), 
%               fill = "red", alpha = .15) +
%   geom_point() + geom_line() + 
%   ylab("Estimated catchability at age") +
%   facet_wrap(~ qname, scales = "free_x")
% @
%
% as before the data, coefficients and variance covariance are all available via
%
% <<results = "hide", eval = FALSE>>=
% coef(qmod)
% vcov(qmod)
% getX(qmod)
% as.data.frame(qmod)
% @
%
% to use them in you need to combine the coefficients together into one `FLPar` (array)
%
% <<eval = FALSE>>=
% b <- do.call(rbind, coef(qmod))
% Sigma <- vcov(qmod)
% X <- getX(qmod)
% dat <- as.data.frame(qmod)
% dat$fit <- exp(c(X %*% b))# + dat$data.centering)
% xyplot(fit ~ age | cname, 
%        subset = year == 2016, data = dat,
%        type = c("l", "p"))
% @
%
% \subsection{Predict}
%
% Predict simply computes the quantities of interest using the estimated coefficients and the design matrix of the model.
%
% <<>>=
% fit.pred <- predict(fit)
% lapply(fit.pred, names)
% @
%
% \subsection{Simulate}
%
% Simulate uses the variance-covariance matrix computed from the Hessian returned by ``ADMB`` and the fitted parameters, to parametrize a multivariate normal distribution. The simulations are carried out using the method `mvrnorm()` provided by the R package [MASS](http://cran.r-project.org/web/packages/MASS/). `r fign('sim')` shows a comparison between the estimated values and the medians of the simulation, while `r fign('sim2')` presents the stock summary of the simulated and fitted data.
%
% <<>>=
% fits <- simulate(fit, 100)
% flqs <- FLQuants(sim=iterMedians(stock.n(fits)), det=stock.n(fit))
% @
%
% <<sim, fig.cap="Median simulations VS fit">>=
% xyplot(data~year|age, groups=qname, data=flqs, type="l",
%        scales=list(y=list(relation="free", draw=FALSE)),
%        auto.key=list(points=FALSE, lines=TRUE, columns=2),
%        par.settings=list(superpose.line=list(col=c("gray35", "black")),
%        strip.background=list(col="gray90")), ylab="")
% @
%
% <<sim2, fig.cap="Stock summary of the simulated and fitted data">>=
% stks <- ple4 + fits
% plot(stks)
% @
%
\section{The statistical catch-at-age stock assessment framework with MCMC}\label{sec:mcmc}

The previous methods were demonstrated using the maximum likelihood estimation method. However, \code{ADMB} can also use MCMC methods to fit the model. This section shows how the \code{sca} methods interface with ADMB to use the MCMC fits. For this section we'll use the hake assessment in mediterranean areas (gsa) 1, 5, 6 and 7.

The likelihood estimate is:


<<>>=
# ll
data(hke1567)
data(hke1567.idx)
fmod <- ~s(age, k = 4) + s(year, k = 8) + s(year, k = 8, by = as.numeric(age == 0)) + s(year, k = 8, by = as.numeric(age == 4))
qmod <- list(~I(1/(1 + exp(-age))))
fit <- sca(hke1567, hke1567.idx, fmodel=fmod, qmodel=qmod)
fit <- simulate(fit, 1000)
@

To run the MCMC method, one needs to configure a set of arguments, which is done by creating a \code{SCAMCMC} object. For details on the MCMC configuration in \code{ADMB} visit the ADMB website.

<<>>=
# mcmc
mc <- SCAMCMC()
# check the default pars
mc
@

Defaults for now are ok, so lets fit the model. Note that the argument \code{fit} must be set to \code{MCMC} and the argument \code{mcmc} takes the \code{SCAMCMC} object. A major check when running MCMC is the acceptance rate, which should be around 0.3. This is a rule of thumb, for more information read the (extensive) literature on MCMC. The slot \code{fitSumm} stores that information.

<<>>=
# fit the model
fitmc00 <- sca(hke1567, hke1567.idx, fmodel=fmod, qmodel=qmod, fit = "MCMC", mcmc=mc)
# check acceptance rate
fitSumm(fitmc00)
@

<<>>=
plot(hke1567 + fitmc00)
@

As mentioned above \code{ADMB} has several options for MCMC. Here we demonstrate one of them, \code{mcprobe} which sets a fat-tailed proposal distribution, as an example of how to use the \code{SCAMCMC} objects.

<<>>=
mc <- SCAMCMC(mcprobe=0.45)
fitmc01 <- sca(hke1567, hke1567.idx, fmodel=fmod, qmodel=qmod, fit = "MCMC", mcmc=mc)
@

All fits together

<<>>=
plot(FLStocks(ll=hke1567 + fit, mc=hke1567 + fitmc00, mc_alt=hke1567 + fitmc01))
@

\subsubsection{Diagnostics with CODA}

We use the package \code{CODA} to run the diagnostics on MCMC fits. One needs to convert the \aFa output into a \code{mcmc} CODA object over which several diagostics can be ran. The mcmc object is a matrix with the parameters (row = iters, cols= pars).

Common diagnostics for MCMC chains is to look at the burn-in period and auto-correlation. The first can be dealt by droping an initial set of iterations, which is done using the function \code{burnin}. The latter is set by the parameter \code{mcsave N}, which defines the iteration's saving rate (the inverse of the thinning rate). This is the rate at which samples of the parameters are saved, such that thinning is effectively discarding draws.

Next fit will run 1000 iterations and save every iter (\code{mcsave=1}). 

<<>>=
library(coda)
mc <- SCAMCMC(mcmc=1000, mcsave=1)
fitmc02 <- sca(hke1567, hke1567.idx, fmodel=fmod, qmodel=qmod, fit = "MCMC", mcmc=mc)
fitmc02.mc <- FLa4a::as.mcmc(fitmc02)
@

The autocorrelation plots will show the very strong correlation across samples, which we want to avoid. ~\ref{fig:acf01} shows autocorrelation for the first parameter.

<<acf01>>=
acf(fitmc02.mc[,1])
@

Ploting the chain for the parameter clearly shows the autocorrelation but also the burnin phase, where there's no information about the parameter. These iterations must to be dropped. 

<<chain01>>=
xyplot(fitmc02.mc[,1])
@

It's also important to check if the distribution of the parameters is normal, which can be done with the \code{densityplot}:

<<dens01>>=
densityplot(fitmc02.mc[,1])
@

Another interesting diagnostic is the Geweke-Brooks Z-score check. This diagnostic indicates if the first and last part of a sample from a Markov chain may not be drawn from the same distribution. It's useful to decide if the first few iterations should e discarded.

<<gew01>>=
geweke.plot(fitmc02.mc[,1])
@

It's clear from the above diagnostics that a burnin phase of about 200 iterations should be considered. With relation to thining one needs to try several values until no autocorrelation exits. 

Next fit will run 10000 iterations and save every iter (\code{mcsave=10}), so that the same 1000 iters are generated by the method. 

<<>>=
mc <- SCAMCMC(mcmc=10000, mcsave=10)
fitmc03 <- sca(hke1567, hke1567.idx, fmodel=fmod, qmodel=qmod, fit = "MCMC", mcmc=mc)
fitmc03.mc <- FLa4a::as.mcmc(fitmc03)
@

The autocorrelation plots still shows a strong correlation across samples, although less than in the previous model.

<<>>=
acf(fitmc03.mc[,1])
@

Next fit will run 100000 iterations and save every iter (\code{mcsave=100}), so that the same 1000 iters are generated by the method. Autocorrelation is much weaker, could still be reduced by increasing \code{mcsave}. 

<<>>=
mc <- SCAMCMC(mcmc=100000, mcsave=100)
fitmc03 <- sca(hke1567, hke1567.idx, fmodel=fmod, qmodel=qmod, fit = "MCMC", mcmc=mc)
fitmc03.mc <- FLa4a::as.mcmc(fitmc03)
@

Next fit will run 200000 iterations and save every iter (\code{mcsave=200}). 

<<>>=
mc <- SCAMCMC(mcmc=200000, mcsave=200)
fitmc03 <- sca(hke1567, hke1567.idx, fmodel=fmod, qmodel=qmod, fit = "MCMC", mcmc=mc)
fitmc03.mc <- FLa4a::as.mcmc(fitmc03)
@

<<>>=
acf(fitmc03.mc[,1])
@

All diagnostics improved with the new thining rate although some other improvements can be done. Note this diagnostics should b echecked for all parameters. For the sake of space the demonstration uses only on the first.

<<>>=
xyplot(fitmc03.mc[,1])
@

<<>>=
densityplot(fitmc03.mc[,1])
@

<<>>=
geweke.plot(fitmc03.mc[,1])
@

The manual "A Guide for Bayesian Analysis in AD Model Builder" by Cole C. Monnahan, Melissa L. Muradian and Peter T. Kuriyam describe and explain a larger group of arguments that can be set when running MCMC with ADMB, which the engine \aFa uses. 

\subsubsection{Confidence interval coverage and MCMC setup}

\subsubsection{Probabilistic assessment (RH code)}

\section{Advanced features}\label{sec:advfeat}

Not sure this section will stay or subsections will be merged in other sections ...

\subsection{Assessing the coverage of confidence intervals}

\subsection{Propagate natural mortality uncertainty}

% In this section we give an example of how uncertainty in natural mortality, set up using the \code{m()} method and the class \code{a4aM} (see chapter XX), is propagated through the stock assessment. We'll start by fitting the default model to the data.
%
% <<>>=
% data(ple4)
% data(ple4.indices)
% fit <- sca(ple4, ple4.indices)
% @
%
% Using \aFa methods we'll model natural mortality using a negative exponential model by age, Jensen's estimator for the level and a constant trend with time. We include multivariate normal uncertainty using the \code{mvrnorm()} method and create 25 iterations.
%
% <<>>=
% nits <- 25
%
% shape <- FLModelSim(model=~exp(-age-0.5))
% level <- FLModelSim(model=~k^0.66*t^0.57, params = FLPar(k=0.4, t=10),
%                     vcov=matrix(c(0.002, 0.01,0.01, 1), ncol=2))
% trend <- FLModelSim(model=~b, params=FLPar(b=0.5), vcov=matrix(0.02))
%
% m4 <- a4aM(shape=shape, level=level, trend=trend)
% m4 <- mvrnorm(nits, m4)
% range(m4)[] <- range(ple4)[]
% range(m4)[c("minmbar","maxmbar")]<-c(1,1)
% flq <- m(m4)[]
% quant(flq) <- "age"
% stk <- propagate(ple4, nits)
% m(stk) <- flq
% @
%
% We fit the same model to the new stock object which has uncertainty in the natural mortality. The assessment is performed for each of the 25 iterations.
%
% <<>>=
% fit1 <- sca(stk, ple4.indices)
% @
%
% And compare the two results (Figure~\ref{fig:mprop}). It's quite easy to run these kind of tests and a large part of our effort is to create the tools to do so.
%
% <<mprop, fig.cap="Stock summary for two M models">>=
% plot(FLStocks("Jensen M with uncertainty"=propagate(ple4, 25)+fit1, "Jensen M"=ple4+fit), key=TRUE)
% @
%
\subsection{Replicating itself (WCSAM) and parallel computing}

% The \href{http://www.ices.dk/news-and-events/symposia/WCSAM-2013}{World Conference on Stock Assessment Methods} promoted a workshop where a large simulation study was used to test the performance of distinct stock assessment models. The first criteria used was that the models should be able to reproduce itself. The process involved fitting the model, simulating observation error using the same model, and refitting the model to each iteration. The final results should be similar to the fitted results before observation error was added (see [Deroba, et.al, 2014](http://icesjms.oxfordjournals.org/content/early/2014/01/18/icesjms.fst237.abstract) for details).
%
% To make things faster one can paralelize the analysis with the \code{parallel} R package. In this example each iteration is a dataset, including surveys, and we'll run one assessment for each iteration. Afterwards the data is pulled back together in an \code{FLStock} object and plotted (Figure~\ref{fig:wcsampar}). Only 10 iterations are run to avoid taking too long. Also note that we're using 4 cores. This parameter depends on the computer being used. These days almost all computers have at least 2 cores.
%
% <<>>=
% # number of iterations
% nits <- 2
% # default fit
% fit0 <- sca(ple4, ple4.indices)
% stk0 <- ple4 + fit0
%
% # simulate and update stock and index
% set.seed(123)
% fit0s <- simulate(fit0, nits, 1234)
% stk0s <- ple4 + fit0s
% idx0s <- ple4.indices
% for(i in 1:length(idx0s)) index(idx0s[[i]]) <- index(fit0s)[[i]]
%
% # paralelize
% library(parallel)
% lst <- mclapply(split(1:nits, 1:nits), function(x){
%   stk <- iter(stk0s, x)
%   idx <- iter(idx0s, x)
% 	out <- try(sca(stk, idx, fit="MP"))
% 	if(is(out, "try-error")) NULL else stk + out
% }, mc.cores=4)
%
% # put things back together
% stks1 <- stk0s
% for(i in 1:nits){
% 	iter(catch.n(stks1), i) <- catch.n(lst[[i]])
% 	iter(stock.n(stks1), i) <- stock.n(lst[[i]])
% 	iter(harvest(stks1), i) <- harvest(lst[[i]])
% }
% catch(stks1) <- computeCatch(stks1)
% stock(stks1) <- computeStock(stks1)
% flsts <- FLStocks(original=stk0, simulation=stk0s, "fit over simulation"=stks1)
% @
%
% <<wcsampar, fig.cap="Replicating the stock assessment model (WCSAM approach) using parallel computing">>= 
% plot(stks3, key=TRUE)
% @
%
\section{A potential stock assessment workflow}\label{sec:saworkflow}





\end{document}