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
```
249

250
## 6. SR model
251
Stock-recruitment model (SR) is the main equation to make forecasting. SR model quality will specify forecast quality.
252

253
### 6.1 Rec/SSB history
254
First you should review rec/ssb relation in whole time series.
255

256
```{r fig.cap="Fig.6.1.1. Log Recruit/SSB numbers (SSB productivity)", echo=FALSE}
257
knitr::include_graphics("output/7_SR_BRP/rec_per_ssb.png")
258
```
259

260
### 6.2 SR relationship
261
Used package: **`r ifelse(config.srmodel.eqsim.use, "EQSIM", "FLSR")`**, model type: **Segreg**
262

263
```{r echo=FALSE, fig.cap="Fig.6.2.1 Stock-recruitment relationship"}
264
if (config.srmodel.eqsim.use) {
265
  knitr::include_graphics("output/7_SR_BRP/eqsim_model.png")
266
} else {
267
    knitr::include_graphics("output/7_SR_BRP/flsr_model.png")
268
}
269
```
270

271
```{r echo=FALSE}
272
if (config.srmodel.eqsim.use) {
273
 sprintf("Excluded years: %s", paste(config.srmodel.eqsim.shift_years, collapse = ", ")) 
274
}
275
```
276

277
Pay your attention to regression inflection point after which one model become a straight line. This point on SSB axis show Blim reference point.  
278
In table above this point equals to **b** coefficient in Segreg model estimates.
279

280
```{r echo=FALSE}
281
if (config.srmodel.eqsim.use) {
282
  result.eqsim.srmodel$sr.det %>% dplyr::mutate_if(is.numeric, funs(as.character(signif(., 5)))) %>%
283
    knitr::kable(., caption = "Table 6.2.1. EQSIM SRR best simulation estimates")
284
} else {
285
  model(result.sr.model)
286
  params(result.sr.model)
287
}
288
```
289

290
## 7. BRP - biological reference points
291
There is 2 way to estimate biological reference point:
292

293
- Simple solution based on rough assumption of Blim reference point (as 20-39% of Bvirgin)
294
- Complex solution based on SR relationship values
295

296
Fmsy estimate can be found only from SRR simulations (in FLSR or EQSIM package).
297

298
Current stock status can be estimated using reference point proportions:
299

300
- Fcurrent/Fmsy will show if stock is in overfishing status (if Fcur/Fmsy > 1)
301
- Fbar/Fmsy will show if stock is in overfishing at latest 3year average 
302
- SSBcurrent/Bpa will show if precautionary point on biomass is reached and stock collapse risk exist (if SSBcur/Bpa < 1) or not (if SSBcur/Bpa > 1)
303
- SSBcurrent/Blim will show if limit reference point on biomass is reached and stock now is collapsing (SSBcurrent/Blim < 1) or not
304

305
### 7.1 Stock reference point
306
In this section current reference point are presented. This reference point obtained from XSA fit results.
307

308
Fbar = `r round(mean(fbar(window(result.xsa.stock, start=lastYear-2))),3)`, F`r lastYear`=`r round(as.numeric(fbar(window(result.xsa.stock, start=lastYear, end=lastYear))),3)`
309

310
SSBbar=`r toString(round(mean(ssb(window(result.xsa.stock, start=lastYear-2))),1))`, SSB`r lastYear`=`r toString(round(as.numeric(ssb(window(result.xsa.stock, start=lastYear, end=lastYear))),1))`
311

312

313
### 7.2 Manual BRP estimation
314
There is a **manual** section of BRP estimates. This approach is very dirty and inaccurate but it is enought to justify BRP for Russian decision makers. 
315

316
Blim estimated in assumption that `max(SSB[2000-2019])` is relevant to Bvirgin, than:
317

318
$$
319
\begin{aligned}
320
B_{lim} = 0.2 * (SSB_{max} + C_{i})
321
\end{aligned}
322
$$
323
Bpa estimation performed using Blim estimate and uncertainty justification in next way:
324
$$
325
\begin{aligned}
326
\sigma = \sqrt{\frac{1}{N-1} \sum_{i=1}^N (F_i - \overline{F})^2}   \; [1] , \\
327
B_{pa} = B_{lim} * e^{+1.645*\sigma} \; [2] 
328
\end{aligned}
329
$$
330
where F - fishing mortality estimate by N year vector.
331

332
Manual estimate results: **Blim**=`r toString(round(result.brp.blim,1))`, **Bpa**=`r toString(round(result.brp.bpa,1))`, sigma = `r toString(round(result.xsa.sigma,2))`
333

334
```{r echo=F, fig.cap="Fig. 7.2.1 XSA SSB estimate and manual BRP"}
335
knitr::include_graphics("output/7_SR_BRP/ssb_blim_bpa_manual.png")
336
```
337

338
```{r echo=F, fig.cap="Fig. 7.2.2 XSA F estimate and manual F0.1"}
339
knitr::include_graphics("output/7_SR_BRP/fbar_f01_manual.png")
340
```
341

342
### 7.3 EQSIM BRP estimation
343
In this section the estimates from EQSIM equilibrium simulation are provided. 
344

345
1. EQSIM_FIT: Blim = `r toString(round(result.eqsim.srmodel$sr.det["b"],1))`, Bpa = `r toString(round(result.eqsim.srmodel$sr.det["b"]*exp(qnorm(0.95)*result.xsa.sigma),1))`
346
2. EQSIM_RUN: Fmsy = `r round(result.eqsim.brp$refs_interval[,1],2)` conf.int95 = [`r round(result.eqsim.brp$refs_interval[,2],2)`,`r round(result.eqsim.brp$refs_interval[,3],2)`].
347

348
```{r echo=F, fig.cap="Fig. 7.2.2 XSA F estimate and manual F0.1"}
349
knitr::include_graphics("output/7_SR_BRP/ssb_blim_bpa_srmodel.png", dpi = 180)
350
```
351

352
### 7.4 Exploitation rate limits
353
Based on old ICES advice E = 0.4 rule is still acceptable. Here is exploitation rate table based on equation `E = F/(F+M)`.
354

355
```{r echo=FALSE}
356
as.data.frame(result.er.table) %>%
357
  dplyr::select(year, data) %>% # MASS package override select() function from dplyr ... shit happens
358
  dplyr::rename(year=year, e=data) %>%
359
  knitr::kable(., digits = 3)
360
```
361

362
```{r echo=FALSE, fig.cap="Fig. 7.4.1. Exploitation rate history"}
363
knitr::include_graphics("output/7_SR_BRP/e_rate04.png")
364
```
365

366
## 8. Short-term forecasts
367
All forecast scenarious sumarized there. Choice your best one. All forecasts based on geomean SRR model (not segmented regression used for BRP estimation).
368

369
Geomean recruit numbers approximation (Geomean SRR) on `r year.terminal+1`-`r year.terminal+3`: `r as.numeric(params(result.sr.model)[1])`.
370

371
```{r echo=F, fig.cap="Fig. 8.1 Forecasts scenarious summary"}
372
knitr::include_graphics("output/8_Forecasts/short_all_in_one.png")
373
```
374

375
```{r echo=F}
376
tmp.f <- window(result.forecast.all, start=year.terminal)
377
nm <- names(tmp.f)
378
y <- (year.terminal):(year.terminal+5)
379
ssb <- list()
380
catch <- list()
381
for (i in 1:length(tmp.f)) {
382
  ssb[[i]] <- as.vector(ssb(tmp.f[[i]]))
383
  catch[[i]] <- as.vector(catch(tmp.f[[i]]))
384
}
385

386
df.ssb <- data.frame(year = y, 
387
                 ssb[[1]],
388
                 ssb[[2]],
389
                 ssb[[3]],
390
                 ssb[[4]],
391
                 ssb[[5]]
392
                 )
393

394
names(df.ssb) <- c("year", nm)
395

396
knitr::kable(df.ssb, caption = "Table 8.1 SSB forecast at different F scenarious")
397

398
df.catch <- data.frame(year = y, 
399
                 catch[[1]],
400
                 catch[[2]],
401
                 catch[[3]],
402
                 catch[[4]],
403
                 catch[[5]]
404
                 )
405
names(df.catch) <- c("year", nm)
406

407
knitr::kable(df.catch, caption = "Table 8.2 Catch forecast at different F scenarious")
408

409
```
410

411

412
## Appendix list
413
Here are raw output from all analysis stage printed in original form.
414

415

416
### Appendix 1. Input data: stock, index, control
417
```{r echo=FALSE}
418
print(data.stock)
419
print(data.indices[[1]])
420
print(data.indices[[2]])
421
print(config.tuning.xsa)
422

423
```
424

425
### Appendix 2. XSA diagnostics
426
```{r echo=FALSE}
427
print(diagnostics(result.xsa.fit))
428
```
429

430
### Appendix 3. XSA and stock result
431
```{r echo=FALSE}
432
print("SSB estimate")
433
print(ssb(result.xsa.stock))
434
print("Recruit numbers")
435
print(rec(result.xsa.stock))
436
print("Stock numbers")
437
print(stock.n(result.xsa.stock))
438
```