07-diagnostics.Rmd

---
---
---

# 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.

## Residuals

Residuals are a ubiquos metrics to check quality of a fit. For `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 `residuals()` method will compute these residuals and generate a object which can be plotted using a set of packed methods. The argument `type` will allow the user to chose which residuals will be computed. By default the method computes standardized residuals.

```{r}
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$.

```{r, 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.

```{r, 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 (`bubble()`) are shown in Figure \@ref(fig:bub). It shows the same information as Figure \@ref(fig:res) but in a multivariate perspective.

```{r, 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.

```{r, 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 `fit='pearson'` (Figure \@ref(fig:resp)).

```{r, 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 `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.

```{r, 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)
```

## 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 `a4a` framework implements a visualization of the fit's predictive skill for both catch-at-age and abundance indices. These are generated by the method `plot()` with the fit object and a `FLStock` (Figure \@ref(fig:selplt)) or `FLIndices` (Figure \@ref(fig:idxplt)) object as arguments.

```{r, selplt, fig.cap="Predict and observed catch-at-age"}
plot(fit, ple4)
```

```{r, idxplt, fig.cap="Predict and observed abundance-at-age"}
plot(fit, ple4.indices)
```

## 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[^07-diagnostics-1]. The implementation of this diagnopstics is done through the method `computeCatchDiagnostics()`, which can be visualized with `plot()` (Figure \@ref(c_d)).

[^07-diagnostics-1]: 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.

```{r, catchdiag, fig.cap="Diagnostics for age aggregated catch in weight"}
c_d <- computeCatchDiagnostics(fit, ple4)
plot(c_d)
```

## Fit summary, information and cross-validation metrics

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

```{r}
fitSumm(fit)
```

Information criteria based metrics are reported with the methods:

```{r}
AIC(fit)
BIC(fit)
```

## The package a4adiags

The package `a4adiags` contains some additional diagnostics based on the `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.

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

```{r, idxrunstest, fig.cap="Runstest for the abundance index"}
plotRunstest(fit, mut09.idx, combine = F) + theme_bw() + facet_wrap(~age)
```

```{r, 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.

## Residuals and submodels misspecifiction

### The "mean" model

To start the analysis we'll fit a "mean" model, where all submodels will be set to an overall average, by using the $\sim 1$ formula. This will be our reference model to see how adding age and year effects will show up in the diagnostic tools, in particular in the residuals.

```{r}
data(hke1567)
data(hke1567.idx)
fit01 <- sca(hke1567, hke1567.idx, fmod=~1, qmod=list(~1), srmod=~1, vmod=list(~1, ~1),  n1mod=~1)
res01 <- residuals(fit01, hke1567, hke1567.idx)
```

The common residuals plot clearly shows a trend across ages (Figure \@ref(fig:meanresbyage)) for both datasets.

```{r, meanresbyyear, fig.cap="Mean fit residuals by year)"}
plot(res01)
```

Which is even clearer when plotting the residuals by age across years.

```{r, meanresbyage, fig.cap="Mean fit residuals by age)",  out.extra = "angle=90", out.height="18cm", out.width="25cm", fig.height=8, fig.width=12}
plot(res01, auxline="l", by="age")
```

### The age effects

These models will introduce age effects in the fishing mortality submodel and catchability submodel. First in the fishinf mortality submodel.

```{r}
fit02 <- sca(hke1567, hke1567.idx, fmod=~factor(age), qmod=list(~1), srmod=~1, vmod=list(~1, ~1),  n1mod=~1)
res02 <- residuals(fit02, hke1567, hke1567.idx)
```

The residuals plot now shows catch at age residuals less stagered, reflecting the modelling of the age effect.

```{r, fageresbyyear, fig.cap="f age effect fit residuals by year)"}
plot(res02)
```

The residuals plot by age shows the same outcome.

```{r, fageresbyage, fig.cap="f age effect fit residuals by age)",  out.extra = "angle=90", out.height="18cm", out.width="25cm", fig.height=8, fig.width=12}
plot(res02, auxline="l", by="age")
```

Follwed by the same addition to the catchability model.

```{r}
fit03 <- sca(hke1567, hke1567.idx, fmod=~1, qmod=list(~factor(age)), srmod=~1, vmod=list(~1, ~1),  n1mod=~1)
res03 <- residuals(fit03, hke1567, hke1567.idx)
```

```{r, qageresbyyear, fig.cap="q age effect fit residuals by year)"}
plot(res03)
```

The residuals plot by age shows the same outcome.

```{r, qageresbyage, fig.cap="q age effect fit residuals by age)",  out.extra = "angle=90", out.height="18cm", out.width="25cm", fig.height=8, fig.width=12}
plot(res03, auxline="l", by="age")
```

Finally both effects are brought together.

```{r}
fit04 <- sca(hke1567, hke1567.idx, fmod=~factor(age), qmod=list(~factor(age)), srmod=~1, vmod=list(~1, ~1),  n1mod=~1)
res04 <- residuals(fit04, hke1567, hke1567.idx)
```

```{r, fqageresbyyear, fig.cap="q age effect fit residuals by year)"}
plot(res04)
```

The residuals plot by age shows the same outcome.

```{r, fqageresbyage, fig.cap="q age effect fit residuals by age)",  out.extra = "angle=90", out.height="18cm", out.width="25cm", fig.height=8, fig.width=12}
plot(res04, auxline="l", by="age")
```

### The fishing mortality year model

This model will introduce an year effect in the fishing mortality submodel on top of the F age effect added before.

```{r}
fit05 <- sca(hke1567, hke1567.idx, fmod=~factor(age) + factor(year), qmod=list(~1), srmod=~1, vmod=list(~1, ~1),  n1mod=~1)
res05 <- residuals(fit05, hke1567, hke1567.idx)
```

The residuals plot now shows catch at age residuals stagered as before. The year trends are less pronounced although, because the data doesn't have a very strong year effect, it's less clear than when modelling the age effect.

```{r, fyearresbyyear, fig.cap="f year effect fit residuals by year)"}
plot(res05)
```

The residuals plot by age shows the same outcome.

```{r, fyearresbyage, fig.cap="f year effect fit residuals by age)",  out.extra = "angle=90", out.height="18cm", out.width="25cm", fig.height=8, fig.width=12}
plot(res05, auxline="l", by="age")
```

We can see now that the residuals show a lot less patterns than before. There's still some issues, the survey catchability seems to have an year trend. However the model is not fully specified yet, stock recruitment is modelled as constant over time, the initial population abundance is also modelled as a constant as well as the variance models.

### The initial year population abundance model, aka N1

This model will introduce an age effect in the population abundance in the first year of the time series. This model sets the n-at-age in the first year of the time series, which is needed due to the lack of previous data to reconstruct those cohorts.

```{r}
fit06 <- sca(hke1567, hke1567.idx, fmod=~factor(age) + factor(year), qmod=list(~factor(age)), srmod=~1, vmod=list(~1, ~1),  n1mod=~factor(age))
res06 <- residuals(fit06, hke1567, hke1567.idx)
```

The residuals plot now shows catch at age residuals stagered as before. The year trends are less pronounced although, because the data doesn't have a very strong year effect, it's less clear than when modelling the age effect.

```{r, n1arresbyyear, fig.cap="f year effect fit residuals by year)"}
plot(res06)
```

The residuals by age (Figure \@ref(fig:n1resbyage)) the residuals' improvement in the first year of the catch at age time series (bottom left plots).

```{r, n1resbyage, fig.cap="f year effect fit residuals by age)",  out.extra = "angle=90", out.height="18cm", out.width="25cm", fig.height=8, fig.width=12}
plot(res06, auxline="l", by="age")
```

### The stock recruitment submodel

In this example we'll simply add a model to allow recruitment to vary over time and we'll see how to track potential improvements in the residuals.

```{r}
fit07 <- sca(hke1567, hke1567.idx, fmod=~factor(age) + factor(year), qmod=list(~factor(age)), srmod=~factor(year), vmod=list(~1, ~1),  n1mod=~factor(age))
res07 <- residuals(fit07, hke1567, hke1567.idx)
```

The residuals plot by year are very useful to see the effect of adding a varying stock recruitment model. The year trends present in previous models are not absent. Recruitment variability when left unmodelled was being picked up by trends in the survey catchability and catch at age. And due to the cohort dynamics underlying the catch at age model, where propagating into other ages' estimates.

```{r, srresbyyear, fig.cap="f year effect fit residuals by year)"}
plot(res07)
```

```{r, srresbyage, fig.cap="f year effect fit residuals by age)",  out.extra = "angle=90", out.height="18cm", out.width="25cm", fig.height=8, fig.width=12}
plot(res07, auxline="l", by="age")
```

### The variance submodel

Finally, we're testing the variance submodel, specifically the catch at age variance model. We won't dig into the catchability variance model though. It's common to accept that a scientific survey following a well designed sampling protocol will have equal variance across ages since no preferential areas are sampled.

```{r}
fit10 <- sca(hke1567, hke1567.idx, fmod=~factor(age) + s(year, k=10), qmod=list(~factor(age)), srmod=~s(year, k=10), vmod=list(~factor(age), ~1),  n1mod=~factor(age))
res10 <- residuals(fit10, hke1567, hke1567.idx)
```

To see what's happening with the variance model one can use predict to plot the different models fitted.

```{r, vagepredbyage, fig.cap="Variance models for catch at age"}
flqs <- FLQuants(mod10=predict(fit10)$vmodel$catch[,"2022"], mod07=predict(fit07)$vmodel$catch[,"2022"])
xyplot(data~age, data=flqs, group=qname, type="l", auto.key=T)
```

To see the effect these models have on the estimated quantities one can look at the variance of the estimates:

```{r, vage, fig.cap="Estimates of population abundance with different variance models"}
flqs <- FLQuants(mod10=catch.n(hke1567+simulate(fit10, nsim=500))[,"2022"], mod07=catch.n(hke1567+simulate(fit07, nsim=500))[,"2022"])
bwplot(data~qname|factor(age),  data=as.data.frame(flqs), scales="free", auto.key=T)
```

and the usual residuals

```{r, vageresbyyear, fig.cap="f year effect fit residuals by year)"}
plot(res10, auxline="r")
```

```{r, vageresbyage, fig.cap="f year effect fit residuals by age)",  out.extra = "angle=90", out.height="18cm", out.width="25cm", fig.height=8, fig.width=12}
plot(res10, auxline="l", by="age")
```