Observations_SAM_SS3_Krill_481.Rmd

---
title: "APPENDIX 1: Integrated approach to modeling krill population dynamics in Western Antarctic Peninsula. Spatial and ecosystem considerations"
subtitle: "Response to observations provided by the Working Group on WP- SAM 2024/26"
output:
  bookdown::pdf_document2:
    number_sections: false
fig_caption: yes
bibliography: SA_krill.bib
csl: apa.csl
link-citations: yes
toc: false
linkcolor: blue
linestretch: 1.3
header-includes:
- \fontsize{12}{16}
- \selectfont
- \usepackage{lscape}
---
\newpage
\tableofcontents
\newpage

```{r setup1, echo=FALSE}
rm(list = ls())
knitr::opts_chunk$set(eval=TRUE,
                      echo = FALSE,
                      message = FALSE,
                      warning = FALSE,
                      fig.align = 'center',
                      fig.pos = "H",
                      dev = 'jpeg',
                      dpi = 300,
                      tidy.opts=list(width.cutoff=80), 
                      tidy=TRUE)
```

```{r message=FALSE, eval=TRUE}
library(here)
library(kableExtra)
library(ggpubr)
library(tibble)
library(readxl)
library(openxlsx)
library(r4ss)
library(flextable)
dir1<-here("s1")
```

\newpage


# OBSERVATION SAM - 23

The main recommendations were extracted from the WG-SAM report @CCAMLRWGSAM2024.

## Observation 1

*2.3 Noting that the authors recalled a previous independent review of an age-based integrated stock assessment for Antarctic krill which encouraged the development of such assessments (Thomson, 2016), the Working Group highlighted that the other reviewer on that panel noted that a length-based model could be considered due to the sparsity of direct age data (de Lestang, 2016). It also noted that this would avoid the approximations required when converting length data to age data. The Working Group further recalled that a similar comment had been made (WG-SAM-2023, paragraph 4.3) in relation to a pilot Casal2 age-based assessment (WG-SAM-2023/25).*

- **Response**: 

We consider a length-based model as an alternative to the age-based approach, especially given the scarcity of direct age data in krill monitoring. Aging techniques for krill have not been available historically, although some recent progress has been made [@Kilada2017]. In this sense, when reading of hard parts is difficult like in krill, catch-at-length data are more plentiful, as the collection of length information is relatively cheap. Length offers insights into the age structure of the population, as there is a correlation between age and length.  So far, all integrated models for krill have been of the Catch-at-Length type [@Kinzey2015a; @Kinzey2019a; @Wang2021]. This means they are approaches that internally model age through growth parameters and the Age-length-key (ALK). To perform this conversion, various methods have been used, but the AKL matrix is perhaps the most important one [@ICCA2003].

In summary, ALK are generated by aging a sub-sample of a population and used to estimate the age distribution from larger length samples. ALKs assume both aged and measured fish are random samples from the same population and should be applied within the same time period to avoid biases. Seasonal or multi-annual applications require careful justification, as using an ALK from a single period can introduce significant bias. Length-stratified sampling is essential, and the process of developing ALKs is labor-intensive, requiring optimal data collection to ensure accuracy  [@ICCA2003]. Formulae exist to estimate the number of age determinations and length measurements is calculated as;


$$
\frac{\text{Number at age for a length group}}{\text{Number of fish aged in that length group}}
$$

The ALK for the time period is raised to the length distribution for that time period:

$$
\text{Raised numbers at age by length group} = \text{Numbers at length} \times \text{Proportion at age for that length}
$$

If the ALK lacks data for certain length groups, data from those groups may be assigned to adjacent ones in the ALK. This process generates an age-length distribution by gear and time period, but caution is needed to avoid significant biases, especially for larger lengths with a broad age range. The numbers at age are summed across length groups, combining variances from ageing and length sampling to produce the age composition for the specified period. Numbers at age for all gears can then be calculated. Numbers at age for all gears can be calculated as:

$$
\sum N_a \times \left(\frac{W_{ct}}{W_{cs}}\right)
$$
where $\sum N_a$ is the sum of sampled numbers at age, W~ct~ is the total commercial catch weight, and W~cs~ is the
sampled commercial catch weight. Variance due to ageing of numbers at age for all gears can be calculated as:
  
$$
\sum \text{Var}_a
$$
where $\sum Var_a$  is the sum of variances due to ageing. Variance due to length sampling of numbers at age can be calculated as:

$$
\sum \text{Var}_l
$$
where $\sum Var_a$  is the sum of variances due to length sampling. The variances should be raised by:

$$
\frac{W_{ct}}{W_{cs}}
$$
where W~ct~ is the total commercial catch weight, and W~cs~ is the sampled commercial catch weight.

The proportions-at-age are then adjusted to find the best fit between the observed size frequency data and that predicted by the proportion-at-age and the ALK. The majority of ALKs are based on growth parameters according to the most suitable formula for the species. In the case of krill, these parameters are derived from the von Bertalanffy growth curve. In @Mardones2023, using `SS3` stock assessment catch-at-lenght model, we use a von Bertalanffy growth relationship with estimated parameters L~inf~, k, and CV. This length-at-age relationship, mediated by survey catchability and selectivity, was used to compare survey biomasses with model estimates of vulnerable biomass (Eq. (A.4)). Weight-at-length was assumed known and calculated using @Maschette2020 parameters.

```{r include=FALSE, message=FALSE}
base.model1 <- SS_output(dir=dir1,
                        covar=T,
                        forecast=T)
```



```{r AKL}
alk_matrix <- base.model1$ALK[,,1] # Ajusta las dimensiones según tu matriz

# Convertir la matriz en un data.frame
alk_df <- as.data.frame(as.table(alk_matrix))

# Renombrar columnas para mayor claridad
colnames(alk_df) <- c("Length", "TrueAge", "Value")

# Asegúrate de que las columnas Length y TrueAge sean numéricas
alk_df$Length <- as.numeric(as.character(alk_df$Length))
alk_df$TrueAge <- as.numeric(as.character(alk_df$TrueAge))

# Asegúrate de que la columna Value sea numérica
alk_df$Value <- as.numeric(alk_df$Value)

# Crear el plot usando ggplot
ggplot(alk_df, aes(x = TrueAge, y = Length, fill = Value)) +
  geom_tile() +
  scale_fill_gradient(low = "white", high = "blue") +
  labs(x = "Age",
       y = "Length",
       fill = "Value") +
  theme_bw()+
  theme(legend.position = "none")
```

One way to avoid dependence on converting variables such as size selectivity to age implemented in `SS3` is through the use of *platoons*. This approach allows for the exploration of survivorship that varies with size against the age. A value of 1 will not create additional groups. Odd-numbered values (e.g., 3, 5) will divide the overall morph into that many groups, resulting in smaller, larger, and average growth groups. Increasing the number of groups slows down the model's performance, so it's recommended not to exceed a value of 5. The fraction of each morph assigned to each group can be either user-defined or approximated using a normal distribution. When multiple groups are designated, an additional input is required to specify the ratio of variability in size-at-age between groups versus within groups. This ratio is used to distribute the total growth variability. The size-at-age for each group is then calculated as a factor (based on the between-group and within-group variability) times the size-at-age of the central morph, which is determined from the growth parameters specific to that Growth Pattern x Sex combination [@methot2020stock].

The age specific rates of fishing mortality can be linked to length specific selectivity. If growth platoons are invoked, then each platoon can have different age specific F due to interaction of length selectivity with their different size at age (Methot, *com pers*)

In `SS3`, this approach was conditioned as follows;

```{r eval=FALSE, echo=TRUE}
#_N_Growth_Patterns (Growth Patterns, Morphs, Bio Patterns, GP are terms used interchangeably in SS)
1  
#_N_platoons_Within_GrowthPattern
3 
 #_Platoon_between/within_stdev_ratio (no read if N_platoons=1)
0.4
#vector_platoon_dist_(-1_in_first_val_gives_normal_approx)
c(0.237, 0.464, 0.237)
```


```{r platoon, out.width='70%', fig.cap="Representation of platoons in Krill length frequencies"}
knitr::include_graphics('Figs/Plattoons.png')
```

in this way, regarding to observation in SAM -2023 (paragraph 4.3);

*4.3 its usefulness in providing an additional approach to assessing krill stock status, it recalled that the development of an integrated stock assessment for krill was considered desirable within three to five years"* but considering this;

*(i) the impact of using an age-based or a length-based implementation could be assessed, noting that a preservation of the source data would be preferred (i.e. the conversion between length composition and age composition should be considered, with conversions from length to age to length avoided)*

Since we do not have an age data monitoring program, we explored a *platoon* alternative to compare its performance and results with the current length-based model currently in WG-SAM-2024/26 [@Mardones2023].
 


\newpage 

## Observation 2

*It encouraged the authors to provide standard model diagnostics, similar to those presented for toothfish assessments, to facilitate understanding of model performance (e.g., WG-SAM-2023, paragraphs 6.33 and 6.34).*

Regarding @WGSAM2114, they provide diagnostic plots for the 2021 assessment model for Antarctic toothfish (*Dissostichus mawsoni*) in the Ross Sea region, following the recommendations of WG-SAM-2015 (SC-CAMLR-XXXIV 2015 Annex 5), and for which it is recommended to follow the diagnostic structure for this analysis. @WGSAM2114 has several way to test the perfomance of stock assessment in toothfish, among them are;  progrssion catch-at-age, likelihood profile, Maximum posterior density (MPD), Autocorrelation plots (ACF), Pearson residuals by sex, Markov Chain Monte Carlo priors and posteriors distributions and fits to ages compositions.

Some of the analyses presented in @WGSAM2114 are also implemented in @Mardones2023.

- 1. Likelihood profile


```{r likepro1, out.width='70%', fig.cap="Log-likelihood profiles for R0 for the various data components included in Antarctic krill"}
knitr::include_graphics('Figs/Profile_s1.png')
```

- 2. Pearson residuals to fishery and surveys

```{r resfleet, out.width='50%', fig.cap="Pearson residuals to fishery. Closed bubbles are positive residuals (observed > expected) and open bubbles are negative residuals (observed < expected). in Antarctic krill"}
knitr::include_graphics('s1/plots/comp_lenfit_residsflt6mkt0.png')
```

```{r ressurv, out.width='50%', fig.cap="Pearson residuals to survey. Closed bubbles are positive residuals (observed > expected) and open bubbles are negative residuals (observed < expected). in Antarctic krill"}
knitr::include_graphics('s1/plots/comp_lenfit_residsflt7mkt0.png')
```
- 3. Numbers by length 

```{r numsex, out.width='50%', fig.cap="Beginning of year expected numbers at length in (max ~ 2890.5 billion) in Antarctic krill"}
knitr::include_graphics('s1/plots/numbers6_len_sex1.png')
```
- 4. Monte Carlo Prior distributions

```{r marcov, out.width='90%', fig.cap="Each parameter distribution plots (plot 1 of 8). Deviation parameters are not included. Plotting range is equal to input limits on parameters in Antarctic krill"}
knitr::include_graphics(c('s1/plots/parameter_distributions_page1.png',
                          's1/plots/parameter_distributions_page2.png',
                          's1/plots/parameter_distributions_page3.png',
                          's1/plots/parameter_distributions_page4.png',
                          's1/plots/parameter_distributions_page5.png',
                          's1/plots/parameter_distributions_page6.png'))
```


- 5. Autocorrelation recruitment


```{r acf, out.width='90%', fig.cap=" Cross correlation in autoregressive analysis for all scenarios in krill population"}
knitr::include_graphics('index_files/figure-html/unnamed-chunk-46-1.jpeg')
```

- 6. Likelihood values

```{r likeall, out.width='90%', fig.cap="Components of total likelihood for all sources and all tested scenarios"}
knitr::include_graphics('Figs/Likelihoodtotal.png')
```


It is important to note that there is a degree of consistency in how the performance of stock assessment models for toothfish and krill is evaluated. For krill, the primary diagnostics conducted adhere to the "best practices" outlined in @Carvalho2021b work. This recipe of steps for testing model performance is illustrated in Figure \@ref(fig:recipediags).


```{r recipediags, out.width='80%', fig.cap="Conceptual process flow chart illustrating a series of interconnected diagnostic tests recommended when developing a base model."}
knitr::include_graphics('Figs/recipediags.png')
```

While we acknowledge the work done by @WGSAM2114 in toothfish stock assessment, we propose a convergence towards these standard practices, regardless of the assessment platform and the structure of the integrated model.

\newpage

## Observation 3.1

*2.4 While noting that this study constituted interesting and important work, the Working Group noted that some issues needed further consideration including the change of trawl designs over the course of the collection of the data used as inputs to the model (e.g. changes in mesh sizes and the presence of fine-mesh codend inserts)...*


This observation raises implications for fishery selectivity. To develop this we configured  a selectivity block in an alternative scenario (`s9`) and compared its performance with the base model (`s1`). This is configured for all strata in 48.1. Tn this case, selectivity blocks are configured in `SS3` as show this template lines in `control` file:

```{r eval=FALSE, echo=TRUE}
1 #_Nblock_Patterns
2 #_blocks_per_pattern
#_begin and end years of blocks
c(1980,2006,2006,2020)
```

The comparison of the main evaluation results with both scenarios is presented in the estimates of the population variables of biomass (Figure \@ref(fig:comparebio)) and recruitment (Figure \@ref(fig:comparerec)). It can be observed that there are no significant differences, at least in the state variables.

```{r comparebio, out.width='80%', fig.cap="Biomass estimates  in s1 and s9 scenarios."}
knitr::include_graphics('Figs/COM1_9compare4_Bratio_uncertainty.png')
```

```{r comparerec, out.width='80%', fig.cap="Biomass estimates  in s1 and s9 scenarios."}
knitr::include_graphics('Figs/COM1_9compare10_recruits_uncertainty.png')
```
## Observation 3.2

*...and the likely invalid assumption of the Peninsula as being a closed system. It further highlighted the need to discuss the development of a standardized data collection plan to support the ongoing revision of the krill fishery management approach.*

Regarding stock unit issue, in the last online workshop of the SCAR Krill Expert Group (SKEG), conclude that The opinion of the SKEG was that the southern boundary of the ACC (sbACC) serves as an important potential pathway for krill flux into Subarea 48.1, from the Bellingshausen Sea to the northern Antarctic Peninsula, and into Subarea 48.2. The Weddell Front and Antarctic Slope Front were also considered important pathways for krill flux into Subarea 48.2, but with focus in 48.1, with a importance of retention and flux (source to destination relation) identified by SKEG should be taken into account when allocation of catch limits, particularly the importance of retention in the Gerlach Strait, Bransfield Strait and Joinville Island strata, and the directions of the important fluxes [@WGEMM2439] (Figure \@ref(fig:ksh)). 


There are two important issues related to the current stock assessment approach presented in WG-SAM 24/27. The first one, while we understand that there are krill flows to and from Subrea 48.1, the challenge of having time series data to implement a model that covers most of Area 48 becomes difficult. However, we note that as monitoring programs provide the necessary data to implement models of this nature, we are open to a more comprehensive approach in the Southern Ocean. Secondly, the current model implemented in `SS3` considered the spatial heterogeneity of the krill stock structure within Subrea 48.1, as illustrated in the Figure \@ref(fig:ksh), however, a model that involves movement between strata has not yet been developed, and is part of the ongoing work presented at WG-SAM.

```{r ksh, out.width='80%', fig.cap="Graphic representation of Krill flux and retention in management strata proposed in the revised krill fishery management strategy within Subarea 48.1, as indicated by the SKEG WG"}
knitr::include_graphics('Figs/KSH.png')
```
\newpage

## Observation 4

*2.5 Dr S. Kasatkina (Russian Federation) noted that data on krill length and biological composition from catches of fishing vessels will be not suitable for such a modelling approach, and recalled that comparisons of krill length composition from catches of fishing trawls and catches of scientific trawls within the same fishing ground revealed significant differences (WG-ASAM-2021/03). Moreover, Dr Kasatkina noted that there were significant differences in the length composition of catches between fishing vessels and these differences are random in nature, which may be due to both the selective properties of commercial trawls, different fishing methods (continuous and traditional fishing) as well as the efficiency of krill sampling by observers at-sea (WG-ASAM-2021/03; WG-EMM-2024/37). Dr Kasatkina noted that this modelling effort required clarity regarding the interaction between the fishery and krill-dependent predators, which requires regular observations to study the spatial overlap of fishing zones and predator foraging zones and could be accompanied by krill distribution patterns (for example, such complex ecosystem observations were provided on RV Atlantida in 2020; SC-CAMLR-42/07). Dr Kasatkina recalled that the spatial and temporal dynamics of groups from the Bellingshausen and Weddell Seas (Fach et al, 2002; Murphy et al, 2004; WG-EMM-2024/43; WG-EMM-2024/39) and noted that without data on krill transport and standardised acoustic surveys of krill, it is impossible to assess the influence of factors such as spatial heterogeneity and life history parameters on key krill population variables.*

We appreciate the identification of potential issues related to changes in trawl designs and the assumption of the Peninsula as a closed system. We will review the trawl data to assess the impact of changes in mesh sizes and fine-mesh codend inserts on the model inputs.  According to the `C1` form from the database provided by the secretary through the SISO program, it is possible to identify (Figure \@ref(fig:trawl)) that there was a change in technique starting in 2006. 


```{r trawl, out.width='70%', fig.cap="Display of length compositions with different trawl technique use in 48.1 subarea."}
knitr::include_graphics('Figs/trawl.png')
```

This consideration was implemented as a new scenario (`s9`) as previously described and implemented in stock assessment of krill.

\newpage

## Observation 5

*2.6 The Working Group encouraged the authors to provide a progression of model implementations from simple to more complex to facilitate understanding and evaluate the evidence for the model assumptions. It further noted the presence of patterns in the residuals shown in the paper (Figures 5 and 6) which warranted further investigation, as well as the need to assess the realism of some parameter values used in the model. The Working Group also indicated that this work would benefit from taking into consideration recent findings by SKEG regarding the krill stock hypothesis (e.g. WG-EMM-2024/39).*

We agree that a progression from simple to more complex model implementations would facilitate understanding and provide a clearer evaluation of the model assumptions. However, and regarding this observation, a visual inspection of Pearson residuals can lead to a mistaken conclusion about the model *goodness-of-fit*. We have integrated a quantitative element to evaluate this type of test. The plotting options are kept mainly to those provided by `r4ss`. With `SSplotRunstest()`, run tests results are illustrated for mean lengths of size composition data from the  `SS3`model. 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 ‘three-sigma limit’ for that series. We plot this in Figure \@ref(fig:runtest).

```{r runtest, out.width='40%', fig.show='hold',  fig.cap="Runs test plot and Joint residual plot for mean lengths from fits length composition data from different fleets in Antarctic krill."}
sspar(mfrow = c(5, 2), plot.cex = 0.8)
knitr::include_graphics(c('Figs/run1.png',
                          'Figs/run2.png',
                          'Figs/run3.png',
                          'Figs/run4.png',
                          'Figs/run5.png',
                          'Figs/run6.png',
                          'Figs/run7.png',
                          'Figs/run8.png',
                          'Figs/run9.png'))
```
To facilitate automated processing, results from several diagnostic tests can also be called as table. Table \@ref(tab:runtable)

```{r runtable, eval=TRUE, echo=FALSE}
data <- data.frame(
  Index = c("FISHERYBS", "FISHERYEI", "FISHERYGS", "FISHERYJOIN", "FISHERYSSIW", 
            "SURVEYBS", "SURVEYEI", "SURVEYGS", "SURVEYJOIN", "PREDATOR"),
  runs.p = c(0.500, 0.145, 0.338, NA, 0.406, 0.189, 0.148, 0.334, 0.500, 0.599),
  test = c("Passed", "Passed", "Passed", "Excluded", "Passed", 
           "Passed", "Passed", "Passed", "Passed", "Passed"),
  sigma3.lo = c(-0.1816665, -0.2319052, -0.1813395, NA, -0.1476043, 
                -0.2452391, -0.2482065, -0.3723597, -0.5749614, -0.3154527),
  sigma3.hi = c(0.1816665, 0.2319052, 0.1813395, NA, 0.1476043, 
                0.2452391, 0.2482065, 0.3723597, 0.5749614, 0.3154527),
  type = rep("len", 10)
)

kable(data, format = "markdown", 
      col.names = c("Index", "runs.p", "test", "sigma3.lo", "sigma3.hi", "type"),
      caption = "Runs tests results for fits differente fleet used in SS3 krill assessment") %>%
  kable_styling(full_width = FALSE)
```

Regarding the reliability of the parameters, it is important to note that the model is conditioned with input parameters gathered from the literature [@Smith2023a; @Maschette2020; @Kinzey2011], as is shown in Table \@ref(tab:parainit). However, the model calculates the most essential parameters through an iterative process, ensuring that their correlation is carefully managed. This calculation process is specifically designed to avoid dependence on the initial conditioning parameters of the model.

\newpage
```{r parainit, eval=TRUE}

# leo archivos para plotear y hacer tablas
start1 <- SS_readstarter(file = file.path(dir1,
                                          "starter.ss"),
                              verbose = FALSE)
# note the data and control file names can vary, so are determined from the 
# starter file.
dat1 <- SS_readdat(file = file.path(dir1, start1$datfile),
                        verbose = FALSE)
# Read in ctl file. Note that the data fileR object is needed so that SS_readctl
# assumes the correct data structure
ctl1 <-  r4ss::SS_readctl(file = file.path(dir1,
                                    start1$ctlfil),
                        verbose = FALSE,
                        use_datlist = TRUE, 
                   datlist = dat1)
fore1 <- r4ss::SS_readforecast(file = file.path(dir1, 
                                                "forecast.ss"),
                              verbose = FALSE)
# can also read in wtatage.ss for an empirical wt at age model using
# r4ss::SS_readwtatage()

parbio<-ctl1$MG_parms[1:10,c(1:4,7)]
 row.names( parbio)<-c("Nat M",
                       "Lmin", 
                       "Lmax",
                       "VonBert K",
                       "CV young",
                       "CV old", 
                       "Wt a", 
                       "Wt b",
                       "L50%", 
                       "Mat slope")

 SRpar<-ctl1$SR_parms[1:5,c(1:4,7)]
 Qpar<-ctl1$Q_parms[1:2,c(1:4,7)]
 Selpar<-ctl1$size_selex_parms[1:22,c(1:4,7)]
 parInit<-as.data.frame(rbind(parbio,SRpar,Qpar,Selpar))

parInit %>%
  kbl(booktabs = TRUE,
      format = "latex",
      position="ht!",
    caption = "Input parameters for the initial SS3 model of krill. Each parameter line contains a minimum value (LO), maximum value (HI), and initial value (INIT). If the phase (PHASE) for the parameter is negative, the parameter is fixed as input") %>%
  kable_paper("hover", 
              full_width = F)%>%
  kable_styling(latex_options = c("striped"),
                full_width = FALSE,
                font_size=9)%>% 
  pack_rows(index = c("Natural Mortality" = 1,
                        "Growth"= 5,
                        "Length-Weigth Relation" = 2,
                        "Maturity"=2,
                        "Stock-Recruit Relation"=5,
                        "Catchability"=2,
                        "Selectivity"=4))
```

\newpage

## Observation 6

Although it was not recorded in the WGSAM-2024 report, one of the questions raised was related to the representation of population variables, such as the harvest rate.  The historical harvest rate of krill, as depicted in the graph, shows a significant decline from the 1970s through the 1990s, followed by a period of stabilization at much lower levels in the 2000s and 2010s. Initially, the exploitation rate was notably high, particularly in the Bransfield stratum, indicating intense harvesting activities during these early decades. However, after the early 1990s, the exploitation rate decreased sharply, suggesting the implementation of more effective management practices or a reduction in fishing pressure. In recent years, the exploitation rate across all strata has remained relatively low and stable, which may reflect sustainable fisheries management efforts. This trend highlights the importance of continued monitoring and adaptive management to ensure the long-term sustainability of krill populations in the face of environmental changes and evolving fishing practices (Figure \@ref(fig:hr)).

```{r hr, out.width='80%', fig.show='hold',  fig.cap="Harvest rate by strata in time series (1980-2020) in Antarctic krill fishey"}
# Convertir la matriz en un formato largo (long format) para ggplot
exploitation_long <- tidyr::pivot_longer(base.model1$exploitation, 
                                         cols = starts_with("FISHERY"), 
                                         names_to = "Fishery", 
                                         values_to = "Exploitation")

# Crear el gráfico
ggplot(exploitation_long, 
       aes(x = Yr, y = Exploitation, color = Fishery)) +
  geom_point() +
  geom_line() +
   geom_line(aes(y = annual_F), color = "black", linetype = "dashed", size = 1) + 
  scale_color_viridis_d()+
  labs(title = "",
       x = "Year",
       y = "Harvest Rat",
       color = "Strata") +
  theme_bw()+
  xlim(1980,2020)

```

And total Harvest rate could be noted in Table \@ref(tab:hrt) and Harvest Rate by Strata regarding outputs of `SS3` in Table \@ref(tab:hrt2)

\newpage

```{r hrt}
# Calcular el acumulado anual de F_std por año
data <- exploitation_long %>%
  dplyr::filter(Yr<2021) %>% 
  group_by(Yr) %>%
  dplyr::summarize(F_std_accumulated = mean(F_std)) %>% 
  dplyr::select(Yr,  F_std_accumulated) %>% 
  dplyr::rename("Harvest Rate" ="F_std_accumulated" ,
                "Year" = "Yr")

data %>%
  kbl(booktabs = TRUE,
      format = "latex",
      position="ht!",
  caption = "Accumulate Harvest Rate by year in whole 48.1 Subarea") %>%
  kable_paper("hover", 
              full_width = F)%>%
  kable_styling(latex_options = c("striped"),
                full_width = FALSE,
                font_size=9)


```

\newpage

```{r hrt2}
data2 <- base.model1$exploitation %>% 
  dplyr::select(1, 7:11) %>% 
  dplyr::rename("Year" = "Yr")
data2 %>% 
  kbl(booktabs = TRUE,
      format = "latex",
      position="ht!",
      digits = 3,
  caption = "Harvest Rate by strata and year") %>%
  kable_paper("hover", 
              full_width = F)%>%
  kable_styling(latex_options = c("striped"),
                full_width = FALSE,
                font_size=9) %>% 
  pack_rows("Forecast", 44, 44)
```

\newpage

# REFERENCES