sprat-xsa-experimental

1
---
2
title: "Black Sea sprat XSA stock assessment (`r as.numeric(data.stock@range['minyear'])`-`r as.numeric(data.stock@range['maxyear'])`)"
3
author: "Author: Piatinskii M., Azov-black sea branch of VNIRO"
4
date: 'Report build date: `r Sys.time()`'
5
output:
6
  html_document:
7
    toc: true
8
    toc_depth: 3
9
    toc_float: true
10
    theme: default
11
---
12

13
```{r echo=FALSE}
14
year.terminal <- as.numeric(data.stock@range["maxyear"])
15
year.start <- as.numeric(data.stock@range["minyear"])
16
```
17

18
## 1. Input data 
19
Here you can review input data summary information. This input data used to perform XSA model and make short-term forecasting.
20

21
### 1.1 Stock (FLStock)
22

23
```{r echo=TRUE}
24
summary(data.stock)
25
```
26

27
### 1.2 Index (FLIndices)
28

29
```{r}
30
summary(data.indices)
31
```
32

33
### 1.3 Parametrization
34

35
```{r echo=FALSE}
36
for (var in ls())
37
  if (startsWith(var, "config.") & var != "config.tuning.xsa") print(sprintf("%s=%s", var, environment()[[var]]))
38
```
39

40
## 2. Preliminary diagnostic
41
In this section you should review all available diagnostic information before start analysis. If diagnostics show big residuals, different slope of log index regression, retrospective miss convergence - you should **stop** analysis with "bad quality" input data assumption. The procedure should include next steps:
42

43
1. Review shrinkage impact factor: how shrinkage Fse change results?
44
2. Review retrospective stability at diffenrent shrinkage Fse levels
45
3. Review Mohn's-Rho retrospective test results
46
4. Review Survey age-vs-age regression lines: all lines should have the same slope (b coef) to be representative by time vector
47

48
### 2.1 Shrinkage impact factor
49
Shrinkage standard error (Fse, or shrinkage s.e.) is a major model covergence factor for model fitting to observed data. By default, Fse = 0.5. 
50

51
```{r fig.cap="Fig.2.1.1. Full XSA estimates at different Fse levels", echo=FALSE}
52
knitr::include_graphics("output/1_FseChoose/summary_all.png")
53
```
54

55
```{r fig.cap="Fig.2.1.2. Last 5 years XSA estimates at different Fse levels", echo=FALSE}
56
knitr::include_graphics("output/1_FseChoose/summary_short.png")
57
```
58

59
```{r echo=FALSE}
60
d <- data.frame(year = (year.start):(year.terminal), 
61
                fse0.5 = as.vector(ssb(result.fse.summary[[1]])),
62
                fse1.0 = as.vector(ssb(result.fse.summary[[2]])), 
63
                fse1.5 = as.vector(ssb(result.fse.summary[[3]])), 
64
                fse2.0 = as.vector(ssb(result.fse.summary[[4]])), 
65
                fse2.5 = as.vector(ssb(result.fse.summary[[5]]))
66
                )
67

68
d2 <- data.frame(year = (year.start):(year.terminal), 
69
                fse0.5 = as.vector(fbar(result.fse.summary[[1]])),
70
                fse1.0 = as.vector(fbar(result.fse.summary[[2]])), 
71
                fse1.5 = as.vector(fbar(result.fse.summary[[3]])), 
72
                fse2.0 = as.vector(fbar(result.fse.summary[[4]])), 
73
                fse2.5 = as.vector(fbar(result.fse.summary[[5]]))
74
                )
75

76
knitr::kable(d, caption = "Table 2.1.1. SSB estimates at different Fse levels")
77
knitr::kable(d2, caption = "Table 2.1.2. Fbar estimates at different Fse levels")
78
```
79

80
### 2.2 Retrospective stability
81
There is classic restrospective analysis with step-by-step reducing input time series length by 1 year. The main target of retrospective analysis is observe model convergence in luck-of-data terms. 
82

83
```{r fig.cap="Fig.2.2.1. Retrospective at Fse = 0.5", echo=FALSE}
84
knitr::include_graphics("output/2_RetroFse/0.5.png")
85
```
86

87
```{r fig.cap="Fig.2.2.2. Retrospective at Fse = 1.0", echo=FALSE}
88
knitr::include_graphics("output/2_RetroFse/1.png")
89
```
90

91
```{r fig.cap="Fig.2.2.3. Retrospective at Fse = 1.5", echo=FALSE}
92
knitr::include_graphics("output/2_RetroFse/1.5.png")
93
```
94

95

96
```{r fig.cap="Fig.2.2.4. Retrospective at Fse = 2.0", echo=FALSE}
97
knitr::include_graphics("output/2_RetroFse/2.png")
98
```
99

100
```{r fig.cap="Fig.2.2.5. Retrospective at Fse = 2.5", echo=FALSE}
101
knitr::include_graphics("output/2_RetroFse/2.5.png")
102
```
103

104
### 2.3. Mohn's-Rho test
105
The basic ICES (2018) procedure to determine model stability over years - retrospective Mohn rho test. Mohn's rho test calculate relative bias for latest 3 years (in default ICES procedure - 5 years) retrospective variations on scale -1 ... +1. Low negative values of *rho* leads to underestimate factor, high positive values - to overestimation (remark: values lowest -0.4 or higher +0.4 shows high variations and low model stability). So procedure is:
106

107
$$
108
\begin{aligned}
109
relbias = (retro - base) / base \\
110
rho = mean(relbias)
111
\end{aligned}
112
$$
113

114
For SSB and F math approach will be
115

116
$$
117
\begin{aligned}
118
\rho_{SSB} = \frac{1}{n} \sum_{i=-5}^{n=0} \frac{(SSB_{i} - \overline{SSB})}{\overline{SSB}} \\
119
\rho_{Fbar} = \frac{1}{n} \sum_{i=-5}^{n=0} \frac{(Fbar_{i} - \overline{Fbar})}{\overline{Fbar}} \\
120
\end{aligned}
121
$$
122

123
where i - year steps from terminal year (terminal-3, terminal-2, ... terminal).
124

125
```{r echo=FALSE}
126
result.retro.diagnostic.ssb %>%
127
  rownames_to_column(., var="year") %>%
128
  mutate_if(is.numeric, ~round(.,1)) %>%
129
  knitr::kable(., caption = "Table 2.3.1. SSB retrospective values", row.names = TRUE)
130
```
131

132
```{r echo=FALSE}
133
result.retro.diagnostic.f %>%
134
  rownames_to_column(., var="year") %>%
135
  mutate_if(is.numeric, ~round(.,3)) %>%
136
  knitr::kable(., caption = "Table 2.3.2. Fbar retrospective values", row.names = TRUE)
137
```
138

139
**Final** Rho estimates: 
140

141
$$\rho_{SSB} = `r round(result.retro.mohnrho.ssb,3)`$$
142
$$\rho_{Fbar} = `r round(result.retro.mohnrho.fbar,3)`$$
143

144

145
### 2.4 Survey index regressions
146
Surveys index regression diagnostics. The same regressions slope show a good quality of survey tuning. Survey vs survey plot (fig. 2.4.3) indicate fish numbers distribution between surveys. High correlation r2 value leads to good survey compatibility.
147

148
```{r fig.cap="Fig.2.4.1. Survey 1 log index in time vector", echo=FALSE}
149
knitr::include_graphics("output/4_SurveysRegressions/survey_1.png")
150
```
151

152
```{r fig.cap="Fig.2.4.2. Survey 2 log index in time vector", echo=FALSE}
153
knitr::include_graphics("output/4_SurveysRegressions/survey_2.png")
154
```
155

156
```{r fig.cap="Fig.2.4.3. Survey vs survey age-to-age comparision", echo=FALSE}
157
knitr::include_graphics("output/4_SurveysRegressions/summary_all.png")
158
```
159

160

161
## 3. Model final tuning (FLXSA.control)
162
There is a final FLXSA.control object with tuning used in XSA fitting.
163

164
```{r echo=FALSE}
165
print(config.tuning.xsa)
166
```
167

168
## 4. XSA results
169
Here the final short result of XSA model fitting. More detailed information can be found in appendix.
170

171
### 4.1. Basic results
172

173
```{r echo=FALSE}
174
df <- data.frame(year = year.start:year.terminal, 
175
                 ssb.tons = as.vector(ssb(result.xsa.stock)),
176
                 rec = as.vector(rec(result.xsa.stock)),
177
                 fbar = round(as.vector(fbar(result.xsa.stock)), 3)
178
                 )
179

180
knitr::kable(df, caption = "Table 4.1. XSA estimate results")
181

182
```
183

184
Survivour numbers at age distribution is a major estimate in XSA analysis. Based on stock numbers and fishing mortality information XSA calculate SSB, Fbar and Recruitment numbers. Fig. 4.1.1 show the stock numbers at age estimation. Fig 4.1.2 show the fishing mortality at age estimates. Fig. 4.1.3 show calculated summary results based on stock numbers: spawning stock biomass(SSB), recruitment numbers (rec) and Fbar estimates.
185

186
```{r fig.cap="Fig.4.1.1. Stock numbers at age by time vector", echo=FALSE}
187
knitr::include_graphics("output/5_XsaEstimates/stock_numbers_at_age.png")
188
```
189

190
```{r fig.cap="Fig.4.1.2. Fishing mortality at age by time vector", echo=FALSE}
191
knitr::include_graphics("output/5_XsaEstimates/f_at_age.png")
192
```
193

194
```{r fig.cap="Fig.4.1.3. Summary XSA results", echo=FALSE}
195
knitr::include_graphics("output/5_XsaEstimates/summary_all.png")
196
```
197

198
### 4.2. Uncertainty variations
199
In this section XSA results presented under uncertainty investigation. The uncertainty has been identified based on FAO approach as st.dev of fishing mortality destribution. The uncertainty is: 
200
$$
201
\begin{aligned}
202
\sigma = \sqrt{\frac{1}{N-1} \sum_{i=1}^N (F_i - \overline{F})^2}   \; [1]  \\
203
\end{aligned}
204
$$
205
Uncertainty justification (sigma) used to determine variations in ssb, rec and F as linear multiplier:
206
$$
207
\begin{aligned}
208
var[SSB]_{i} = SSB_{i} *  \sigma  \; [2]  \\
209
var[REC]_{i} = REC_{i} *  \sigma  \; [3]  \\
210
var[F]_{i} = F_{i} *  \sigma  \; [4]  \\
211
\end{aligned}
212
$$
213
Then variation used to create norman distribution of estimates based on mean and sd. The final results shown at fig. 4.2.1 
214

215
```{r fig.cap="Fig. 4.2.1. Uncerteinty factor for Rec, SSB, F", echo = FALSE}
216
knitr::include_graphics("output/5_XsaEstimates/summary_uncertainty.png")
217
```
218

219
In plot technique the classic quantile method used for 0.05(lower), 0.5 (mean) and 0.95 (upper) levels.
220

221
## 5. XSA diagnostics
222
In process of XSA model fit iterations survey information used to tune stock data. Observed vs fitted values (with survey) should be reviewed as fitting quality factor. High residual values leads to bias in model. In fact values in range -2.25 < log(residual) < 2.25 are not significant. Values out of this range idicate low survey quality (by single year or age group). Full negative (year or age group) residual indicate high posibility of bias. 
223

224
```{r fig.cap="Fig.5.1. Log residual diagnostics for survey1 vs model fit", echo=FALSE}
225
knitr::include_graphics("output/6_Residuals/summary_survey1.png")
226
```
227

228
```{r fig.cap="Fig.5.2. Log residual diagnostics for survey2 vs model fit", echo=FALSE}
229
knitr::include_graphics("output/6_Residuals/summary_survey2.png")
230
```
231

232
```{r fig.cap="Fig.5.3. Log catchability residuals at age diagnostics", echo=FALSE}
233
knitr::include_graphics("output/6_Residuals/catch_at_age_summary.png")
234
```
235

236
Log-residual values raw output:
237
```{r echo=FALSE}
238
print(result.xsa.fit@index.res)
239
```
240

241
Min/max residuals:
242
```{r echo=FALSE}
243
for (sur in names(result.xsa.fit@index.res)) {
244
  print(paste("Survey: ", sur))
245
  cat(paste("Min log.res =", min(result.xsa.fit@index.res[[sur]]), "\n"))
246
  cat(paste("Max log.res =", max(result.xsa.fit@index.res[[sur]]), "\n"))
247
}
248
```