index.Rmd

---
title: "Integrated approach to modeling to krill population dynamics in Wester Antarctic Peninsula: spatial and ecosystem considerations"
subtitle: "Working Paper to be submitted in a CCAMLR EMM-WG 2024"
author: "Mardones, M; Cárdenas, C., Krüger, L., Santa Cruz, F."
date:  "`r format(Sys.time(), '%d %B, %Y')`"
bibliography: SA_krill.bib
csl: apa.csl
link-citations: yes
linkcolor: blue
output:
  html_document:
    keep_md: true
    toc: true
    toc_deep: 3
    toc_float:
      collapsed: false
      smooth_scroll: false
    theme: cosmo
    fontsize: 0.9em
    linestretch: 1.7
    html-math-method: katex
    self-contained: true
    code-tools: true
editor_options: 
  markdown: 
    wrap: 72
---

```{r setup1, echo=FALSE}
#rm(list = ls())
#set.seed(999)
knitr::opts_chunk$set(echo = FALSE,
                      message = FALSE,
                      warning = FALSE,
                      fig.align = 'center',
                      fig.pos = "H",
                      dev = 'jpeg',
                      dpi = 300,
                      tidy.opts=list(width.cutoff=50), 
                      tidy=TRUE)
#XQuartz is a mess, put this in your onload to default to cairo instead
options(bitmapType = "cairo") 
# (https://github.com/tidyverse/ggplot2/issues/2655)
# Lo mapas se hacen mas rapido
```

```{r message=FALSE, eval=TRUE}
# install.packages("devtools")
# devtools::install_github("r4ss/r4ss", ref="development")
# install.packages("caTools")
# library("caTools")
# install.packages("r4ss")
library(r4ss)
library(here)
#remotes::install_github("PIFSCstockassessments/ss3diags")
library(ss3diags)
library(kableExtra)
library(doParallel) # facilita la ejecución paralela en R
detectCores()
registerDoParallel(8)
library(ggpubr)
library(tibble)
library(openxlsx)
library(ggthemes)
library(forecast)
library(tidyr)
library(mixR)
library(readxl)
```

```{r message=FALSE, warning=FALSE, echo=TRUE}
dir01 <- here("s01") # agreggate data (no spatial diferences)
dir1<-here("s1") # Data strata flishery
dir2<-here("s2") # Same 9 with areas (SubStrata) as fleet. Dif size comoposition and dif CPUE and dif survey length and biomass data by strata
dir3<-here("s3") # whitout S-R
dir4<-here("s4") # 
dir5<-here("s5") # 
dir6<-here("s6") # 
dir7<-here("s7") # 2 set parametres EMM-2024/23 (Mardones)
dir8<-here("s8") # s1 platoons
dir9<-here("s9") # s1 w/ blocks
dirtest<-here("test")
Figs <- here("Figs")# S
```

# OVERVIEW

In a simple way, the core of Stock Synthesis is its population dynamics
model, which represents the dynamics of krill populations over time.
This model incorporates key biological, environmental and predator data
sources. The model is typically formulated using mathematical equations
that describe how these parameters interact to determine the abundance
and distribution of krill in the study area.

## Statistical model (SS3)

Stock Synthesis v.3.30.21 is a widely used software tool for assessing
fish and invertebrate populations, including krill (Euphausia superba)
in the Antarctic Peninsula region. The methodology employed by Stock
Synthesis involves a comprehensive and integrated approach, utilizing
various data sources and modeling techniques to estimate the main
population variables of krill in WAP.

The stock assesment model was configured using Stock Synthesis (SS3
hereafter)[SS3](https://vlab.noaa.gov/web/stock-synthesis) [@Methot2013;
@methot2020stock] with the most updated version (V.3.30.21). SS3 is a
structured age and size stock evaluation model, in the class of models
called *"Integrated stock evaluation analysis model"*. SS3 has a stock
population sub-model that simulates growth, maturity, fecundity,
recruitment, movement, and mortality processes, and observation
sub-models and expected values for different types of data. The model is
coded in `C++` with estimation parameters enabled by automatic
differentiation (ADMB) [@Fournier2012; @Methot2013]. The analysis of
results and outputs uses R tools and the graphical interface of the
*r4ss* and *ss3diags* library [@Taylor2019; @Winker2023].

By integrating data from multiple sources and considering spatial
heterogeneity, the assessment methodology using Stock Synthesis
v.3.30.21 provides a robust framework for evaluating the status of krill
populations in the Antarctic Peninsula region. This information is
essential for supporting management decisions aimed at ensuring the
sustainable use of krill resources in this ecologically sensitive area.

## Parametres

read files

```{r}
# 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()
```

```{r}
parbio<-ctl1$MG_parms[1:10,c(1:3,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:3,7)]
 Qpar<-ctl1$Q_parms[1:2,c(1:3,7)]
 Selpar<-ctl1$size_selex_parms[1:22,c(1:3,7)]
 parInit<-as.data.frame(rbind(parbio,SRpar,Qpar,Selpar))
```


```{r eval=FALSE}
wb <- createWorkbook()

addWorksheet(wb, "parameters")
writeData(wb, "parameters", parInit)

# Guardar el workbook
saveWorkbook(wb, "DataKrill.xlsx", overwrite = TRUE)

```

```{r}
parInit %>%
  kbl(booktabs = T,
      format = "html",
      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",
                                  "condensed"),
                full_width = FALSE,
                font_size=10,
                html_font = "arial")%>% 
  pack_rows(index = c("Mortalidad natural" = 1,
                        "Crecimiento"= 5,
                        "Relación longitud-peso" = 2,
                        "Ojiva de madurez"=2,
                        "Relación stock-recluta"=5,
                        "Capturabilidad"=2,
                        "Selectividad"=4))

```

## Scenarios

In Table 1 we have ten scenarios to test different option in modeling
about main consideration in assessment of krill population.

| Scenario | Description                                                                          |
|:------------:|:---------------------------------------------------------|
|   s01    | Fishery and Survey (AMLR) data, Predator, Environmental aggregate data in 48.1       |
|    s1    | Fishery and Survey (AMLR) data Length, Index, Catch by strata. Predator and Env data |
|    s2    | "s1" without S-R relation                                                            |
|    s3    | "s1" BH S-R relation weak (0.9 steepness)                                            |
|    s4    | "s1" BH S-R relation strong (0.6 steepness)                                          |
|    s5    | "s1" BH S-R relation mid-strong estimated                                            |
|    s6    | "s1" Ricker S-R relation estimated                                                   |
|    s7    | "s1" w/ set of parameters estimated in @EMM-204/32                                   |
|    s8    | "s1" test data weighting                                                             |
|    s9    | "s1" w/ selectivity blocks                                                             |

## Run Models

```{r eval=FALSE, message=F, include=FALSE, echo=TRUE}
# Lista de directorios para correro tordos juntos
directorios <- c("s01", 
                 "s1", 
                 "s2",
                 "s3",
                 "s4",
                 "s5",
                 "s6",
                 "s7",
                 "s8",
                 "s9")  # Agrega aquí todos los nombres de las carpetas que deseas procesar

# Bucle para ejecutar el código en cada directorio
for (dir in directorios) {
  r4ss::run(
    dir = dir,
    exe = "ss_osx",
    skipfinished = FALSE,
    show_in_console = TRUE
  )
}
```

```{r eval=FALSE, message=F, include=FALSE}
#OR with rby separate
#shell(cmd="ss") # run SS to windows
# or 
r4ss::run(
  dir = dir8,
  exe = "ss_osx",
  skipfinished = FALSE, # TRUE will skip running if Report.sso present
  show_in_console = TRUE # change to true to watch the output go past
)
```

Read outputs

```{r message=F, include=FALSE}
#s01
base.model01 <- SS_output(dir=dir01,
                        covar=T,
                        forecast=T)
#s1
base.model1 <- SS_output(dir=dir1,
                        covar=T,
                        forecast=T)
#s2
base.model2 <- SS_output(dir=dir2,
                        covar=T,
                        forecast=T)
#s3
base.model3 <- SS_output(dir=dir3,
                        covar=T,
                        forecast=T)
#s4
base.model4 <- SS_output(dir=dir4,
                        covar=T,
                        forecast=T)
#s5
base.model5 <- SS_output(dir=dir5,
                        covar=T,
                        forecast=T)
#s6
base.model6 <- SS_output(dir=dir6,
                        covar=T,
                        forecast=T)

#s7
base.model7 <- SS_output(dir=dir7,
                        covar=T,
                        forecast=T)

#s8
base.model8 <- SS_output(dir=dir8,
                        covar=T,
                        forecast=T)

#s9
base.model9 <- SS_output(dir=dir9,
                        covar=T,
                        forecast=T)

#stest
base.modeltest <- SS_output(dir=dirtest,
                        covar=T,
                        forecast=T)
```


Just some time series plots

```{r}
SSplotTimeseries(base.modeltest,
                 subplot = 1)
SSplotTimeseries(base.modeltest,
                 subplot = 9)
```
and fits

```{r}
SSplotComps(base.modeltest,
            subplots = 1,
            maxrows = 5,
            maxcols = 5)
```

```{r eval= FALSE, message=F, include=FALSE}
# html
SS_plots(base.modeltest,
         uncertainty=T,
         datplot = T, 
         png=T, 
         aalresids = F,
         btarg=0.75,
         minbthresh=0.20, 
         forecast=T)

```
# RESULTS

### Main Variables poulation

Total biomass

```{r eval=FALSE}
total01 <- as.data.frame(base.model01$timeseries[1:50,c(2,5)]) %>%
  mutate(Serie = "No Spatial")
total1 <- as.data.frame(base.model1$timeseries[,c(2,5)]) %>% 
  mutate(Serie = "Spatial implicit")
total2 <- as.data.frame(base.model2$timeseries[,c(2,5)]) %>% 
  mutate(Serie = "Spatial w/ without S-R relation")
total3 <- as.data.frame(base.model3$timeseries[,c(2,5)]) %>% 
  mutate(Serie = "Spatial w/ S-R relation weak (0.9 steepness)")
total4 <- as.data.frame(base.model4$timeseries[,c(2,5)]) %>% 
  mutate(Serie = "Spatial w/ S-R relation strong (0.6 steepness)")
total5 <- as.data.frame(base.model5$timeseries[,c(2,5)]) %>% 
  mutate(Serie = "Spatial w/ S-R relation mid-strong estimated ")
total6 <- as.data.frame(base.model6$timeseries[,c(2,5)]) %>% 
  mutate(Serie = "Spatial w/ Ricker S-R relation estimated")
total7 <- as.data.frame(base.model7$timeseries[,c(2,5)]) %>% 
  mutate(Serie = "Spatial w/ New set parameters Linf, k and M")


totalbio <- rbind(total01,
              total1,
              total2,
              total3,
              total4,
              total5,
              total6,
              total7)


ggplot(totalbio, aes(x = Yr, y = Bio_all, color=Serie)) +
  geom_point(size=0.7,
             col="black") +
  geom_smooth(method = "loess",
              span=0.75, 
              se=TRUE)+
  labs(title = "Compare variables between scenarios",
       x = "",
       y = "Total biomass (mt)") +
  theme_few() +
  scale_color_viridis_d(option="H",
                        name="Scenarios")+
  theme(axis.text.x = element_text(angle = 90,
                                   hjust = 1))
```

```{r eval=FALSE, message=F, include=FALSE}
SS_plots(base.model7, 
         uncertainty = TRUE,
         datplot = T, 
         png=T, 
         aalresids = F,
         btarg=0.75, 
         minbthresh=0.2, 
         forecast=T)
```

Data used en both (spatial and No spatial models)

```{r }
SSplotData(base.model01, subplot = 1, 
           fleetnames = c("Fishery", "Survey", "Predator"),
           fleetcol = c("blue","green", "red"))

SSplotData(base.model1, 
           subplot = 1,
           pheight = 15)
```


Respecto a los valores y parametros biologicos modelados, los siguientes
graficos identifican los estimadores puntuales del recurso

```{r fig.height=6}
SSplotSelex(base.model01,
            subplots = 1)
SSplotSelex(base.model1,
            subplots = 1)
SSplotSelex(base.model2,
            subplots = 1)
SSplotSelex(base.model3,
            subplots = 1)
SSplotSelex(base.model4,
            subplots = 1)
SSplotSelex(base.model5,
            subplots = 1)
SSplotSelex(base.model6,
            subplots = 1)
SSplotSelex(base.modeltest,
            subplots = 1)
```



```{r eval=FALSE}
SSplotComps(base.model01,
            subplots = 1,
            maxrows = 5,
            maxcols = 5)
SSplotComps(base.model1,
            subplots = 1,
            maxrows = 5,
            maxcols = 5)
SSplotComps(base.model2,
            subplots = 1,
            maxrows = 5,
            maxcols = 5)
SSplotComps(base.model3,
            subplots = 1,
            maxrows = 5,
            maxcols = 5)
SSplotComps(base.model4,
            subplots = 1,
            maxrows = 5,
            maxcols = 5)
SSplotComps(base.model5,
            subplots = 1,
            maxrows = 5,
            maxcols = 5)
SSplotComps(base.model6,
            subplots = 1,
            maxrows = 5,
            maxcols = 5)
SSplotComps(base.model7,
            subplots = 1,
            maxrows = 5,
            maxcols = 5)
```
Heatmap

```{r}
df <- base.modeltest$lendbase

# Filtra el dataframe
df1 <- df %>% 
  dplyr::filter(Pearson<5) %>% 
  dplyr::select(c(1, 6, 19, 16)) %>% 
  dplyr::mutate(Fleet = dplyr::case_when(
    Fleet == 1 ~ "FISHERYBS",
    Fleet == 2 ~ "FISHERYEI",
    Fleet == 3 ~ "FISHERYGS",
    Fleet == 4 ~ "FISHERYJOIN",
    Fleet == 5 ~ "FISHERYSSIW",
    Fleet == 6 ~ "SURVEYBS",
    Fleet == 7 ~ "SURVEYEI",
    Fleet == 8 ~ "SURVEYGS",
    Fleet == 9 ~ "SURVEYJOIN",
    Fleet == 10 ~ "SURVEYSSIW",
    Fleet == 11 ~ "PREDATOR",
    TRUE ~ NA_character_  # Para manejar cualquier valor que no coincida
  ))

# Crea el heatmap
ggplot(df1, aes(x = Yr, y = Bin, fill = Pearson)) +
  geom_tile() +
  scale_fill_gradient2(low = "blue", 
                       mid = "white", 
                       high = "red",
                       midpoint = 0) +
  facet_wrap(Fleet~., ncol=5) + 
  labs(title = "",
       x = "",
       y = "Length",
       fill = "Pearson") +
  theme_bw()+
  theme(axis.text.x = element_text(angle = 90,
                                   hjust = 1))

```
hexagon


```{r eval=FALSE}
# Assuming your data is in a data frame called 'lendbase'
# Filter and select relevant columns
data <- base.model1$lendbase %>%
  dplyr::filter(Pearson<5) %>% 
  dplyr::select(c(1, 6, 19, 16)) %>% 
  dplyr::mutate(Fleet = dplyr::case_when(
    Fleet == 1 ~ "FISHERYBS",
    Fleet == 2 ~ "FISHERYEI",
    Fleet == 3 ~ "FISHERYGS",
    Fleet == 4 ~ "FISHERYJOIN",
    Fleet == 5 ~ "FISHERYSSIW",
    Fleet == 6 ~ "SURVEYBS",
    Fleet == 7 ~ "SURVEYEI",
    Fleet == 8 ~ "SURVEYGS",
    Fleet == 9 ~ "SURVEYJOIN",
    Fleet == 10 ~ "SURVEYSSIW",
    Fleet == 11 ~ "PREDATOR",
    TRUE ~ NA_character_  # Para manejar cualquier valor que no coincida
  )) %>% 
  dplyr::select(Yr, Bin, Pearson) 

# Create the hexagonal heatmap
ggplot(data, aes(x = Yr, y = Bin)) +
  geom_hex(bins=50) +  # Hexagonal bins
  scale_fill_gradient2(low = "blue", 
                       mid = "white", 
                       high = "red", 
                       midpoint = 0) +  # Color scale
  theme_minimal() +  # Minimal theme for cleaner look
  labs(title = "Hexagonal Heatmap of Pearson Residuals",
       x = "Year",
       y = "Bin",
       fill = "Pearson Residual")  # Labels for plot and axes


ggplot(data, aes(x = Yr, y = Bin)) +
  geom_hex(bins = 50) +
  scale_fill_gradient2(low = "blue", 
                       mid = "white", 
                       high = "red", midpoint = median(data$Bin)) +
  theme_minimal() +
  labs(title = "Hexagonal Heatmap of Pearson Residuals",
       x = "Year", y = "Bin", fill = "Pearson Residual")
```

## Diagnosis Base Model

Step to do a good practice in model diagnosis is;

-     1.  Convergence. Final convergence criteria is 1.0e-04

-     2.  Residual (visual and metrics)

-     3.  Retrospective analysis (Mhon Parameter)

-    4.  Likelihood profile

all this framework try to follow recommendations of @Carvalho2021b

## Residual consistency 

```{r}
SSplotRunstest(base.modeltest,
               subplots = "len",
               add=T,
               plot = TRUE,
               plotdir = Figs)
```


```{r eval=FALSE, fig.width=4, fig.height=5}
SSplotRunstest(base.model01,
               subplots = "len",
               add=T)
SSplotRunstest(base.model1,
               subplots = "len",
               add=T)
SSplotRunstest(base.model2,
               subplots = "len",
               add=T)
SSplotRunstest(base.model3,
               subplots = "len",
               add=T)
SSplotRunstest(base.model4,
               subplots = "len",
               add=T)
SSplotRunstest(base.model5,
               subplots = "len",
               add=T)
SSplotRunstest(base.model6,
               subplots = "len",
               add=T)
SSplotRunstest(base.modeltest,
               subplots = "len",
               add=T)
```

```{r}
SSplotRunstest(base.modeltest,
               subplots = "cpue",
               add=T)
```

```{r eval=FALSE}
SSplotRunstest(base.model01,
               subplots = "cpue",
               add=T)
SSplotRunstest(base.model1,
               subplots = "cpue",
               add=T)
SSplotRunstest(base.model2,
               subplots = "cpue",
               add=T)
SSplotRunstest(base.model3,
               subplots = "cpue",
               add=T)
SSplotRunstest(base.model4,
               subplots = "cpue",
               add=T)
SSplotRunstest(base.model5,
               subplots = "cpue",
               add=T)
SSplotRunstest(base.model6,
               subplots = "cpue",
               add=T)
SSplotRunstest(base.modeltest,
               subplots = "cpue",
               add=T)
```


```{r eval=FALSE}
SSplotJABBAres(base.model01,
               subplots = "len",
               add=T)
SSplotJABBAres(base.model1,
               subplots = "len",
               add=T)
SSplotJABBAres(base.model2,
               subplots = "len",
               add=T)
SSplotJABBAres(base.model3,
               subplots = "len",
               add=T)
SSplotJABBAres(base.model4,
               subplots = "len",
               add=T)
SSplotJABBAres(base.model5,
               subplots = "len",
               add=T)
SSplotJABBAres(base.model6,
               subplots = "len",
               add=T)
SSplotJABBAres(base.modeltest,
               subplots = "len",
               add=T)
```

```{r}
SSplotJABBAres(base.model01,
               subplots = "cpue",
               add=T)
SSplotJABBAres(base.model1,
               subplots = "cpue",
               add=T)
SSplotJABBAres(base.model2,
               subplots = "cpue",
               add=T)
SSplotJABBAres(base.model3,
               subplots = "cpue",
               add=T)
SSplotJABBAres(base.model4,
               subplots = "cpue",
               add=T)
SSplotJABBAres(base.model5,
               subplots = "cpue",
               add=T)
SSplotJABBAres(base.model6,
               subplots = "cpue",
               add=T)
SSplotJABBAres(base.modeltest,
               subplots = "cpue",
               add=T)
```


## Relationship Stock-Recruit

```{r eval=FALSE}
l <- base.model1$recruit$SpawnBio
m <- base.model1$recruit$pred_recr
ggplot()+
  aes(l,m)+
  geom_point()+
  geom_smooth(method = "loess", 
              se =F)+
  theme_few()+
  ylim(0, max(base.model1$recruit$pred_recr))
```

## Retrospective analysis

Los análisis retrospectivo, dan cuenta de diferencias de estimación
(sub - sobreestimación) en los patrones entre modelos evaluados.

```{r eval=FALSE}
#one by one
retro(
    dir = dirtest,
    oldsubdir = "",
    newsubdir = "Retrospective",
    years = 0:-4,
    exe = "ss_osx",
    extras = "-nox",
    skipfinished = FALSE)
```


```{r eval=FALSE}
for (dir in directorios) {
  retro(
    dir = dir,
    oldsubdir = "",
    newsubdir = "Retrospective",
    years = 0:-5,
    exe = "ss_osx",
    extras = "-nox",
    skipfinished = FALSE
  )
}
```

```{r echo=FALSE, message=FALSE, warning=FALSE}
#s01
retroModels01 <- SSgetoutput(dirvec=file.path(dir01,
                                            "Retrospective",
                                            paste("retro",0:-5,
                                                  sep="")))

retroSummary01 <- SSsummarize(retroModels01)
# s1
retroModels1 <- SSgetoutput(dirvec=file.path(dir1,
                                            "Retrospective",
                                            paste("retro",0:-5,
                                                  sep="")))

retroSummary1 <- SSsummarize(retroModels1)
#s2
retroModels2 <- SSgetoutput(dirvec=file.path(dir2,
                                            "Retrospective",
                                            paste("retro",0:-5,
                                                  sep="")))

retroSummary2 <- SSsummarize(retroModels2)
#s3
retroModels3 <- SSgetoutput(dirvec=file.path(dir3,
                                            "Retrospective",
                                            paste("retro",0:-5,
                                                  sep="")))

retroSummary3 <- SSsummarize(retroModels3)
#s4
retroModels4 <- SSgetoutput(dirvec=file.path(dir4,
                                            "Retrospective",
                                            paste("retro",0:-5,
                                                  sep="")))

retroSummary4 <- SSsummarize(retroModels4)
#s5
retroModels5 <- SSgetoutput(dirvec=file.path(dir5,
                                            "Retrospective",
                                            paste("retro",0:-5,
                                                  sep="")))

retroSummary5 <- SSsummarize(retroModels5)
#s6
retroModels6 <- SSgetoutput(dirvec=file.path(dir6,
                                            "Retrospective",
                                            paste("retro",0:-5,
                                                  sep="")))

retroSummary6 <- SSsummarize(retroModels6)
#s7
retroModels7 <- SSgetoutput(dirvec=file.path(dir7,
                                            "Retrospective",
                                            paste("retro",0:-5,
                                                  sep="")))

retroSummary7 <- SSsummarize(retroModels7)
```


```{r eval=FALSE, message=FALSE, warning=FALSE}
#stest
retroModelstest <- SSgetoutput(dirvec=file.path(dirtest,
                                            "Retrospective",
                                            paste("retro",0:-3,
                                                  sep="")))

retroSummarytest <- SSsummarize(retroModelstest)


```

```{r echo=FALSE}
#save(retroSummary, retroModels, file="retro5.Rdata")

retro01 <- SSplotRetro(retroSummary01,
            add=T,
            forecast = F,
            legend = T,
            verbose=F)
retro1 <- SSplotRetro(retroSummary1,
            add=T,
            forecast = F,
            legend = T,
            verbose=F)
retro2 <- SSplotRetro(retroSummary2,
            add=T,
            forecast = F,
            legend = T,
            verbose=F)
retro3 <- SSplotRetro(retroSummary3,
            add=T,
            forecast = F,
            legend = T,
            verbose=F)
retro4 <- SSplotRetro(retroSummary4,
            add=T,
            forecast = F,
            legend = T,
            verbose=F)
retro5 <- SSplotRetro(retroSummary5,
            add=T,
            forecast = F,
            legend = T,
            verbose=F)

retro6 <- SSplotRetro(retroSummary6,
            add=T,
            forecast = F,
            legend = T,
            verbose=F)


retro7 <- SSplotRetro(retroSummary7,
            add=T,
            forecast = F,
            legend = T,
            verbose=F)

retrotest <- SSplotRetro(retroSummarytest,
            add=T,
            forecast = F,
            legend = T,
            verbose=F)

```


```{r eval=FALSE}
retrotest <- SSplotRetro(retroSummarytest,
            add=T,
            forecast = F,
            legend = T,
            verbose=F)
```


## Hindcast Cross-Validation and prediction skill

Implementing the Hindcast Cross-Validation (HCxval) diagnostic in Stock
Synthesis requires the same model outputs generated by
`r4ss:SS_doRetro()`. As a robust measure of prediction skill, we
implemented the mean absolute scaled error (MASE). In brief, the MASE
score scales the mean absolute. Regarding (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.

```{r}
hci = SSplotHCxval(retroSummary1, 
                   add = T, 
                   verbose = F, 
                   legendcex = 0.7)
```


```{r eval=FALSE}
hci = SSplotHCxval(retroSummary01, 
                   add = T, 
                   verbose = F, 
                   legendcex = 0.7)
hci = SSplotHCxval(retroSummary1, 
                   add = T, 
                   verbose = F, 
                   legendcex = 0.7)
hci = SSplotHCxval(retroSummary2, 
                   add = T, 
                   verbose = F, 
                   legendcex = 0.7)
hci = SSplotHCxval(retroSummary3, 
                   add = T, 
                   verbose = F, 
                   legendcex = 0.7) 
hci = SSplotHCxval(retroSummary4, 
                   add = T, 
                   verbose = F, 
                   legendcex = 0.7)
hci = SSplotHCxval(retroSummary5, 
                   add = T, 
                   verbose = F, 
                   legendcex = 0.7)
hci = SSplotHCxval(retroSummary6, 
                   add = T, 
                   verbose = F, 
                   legendcex = 0.7)
hci = SSplotHCxval(retroSummary7, 
                   add = T, 
                   verbose = F, 
                   legendcex = 0.7)
```

## Kobe (status)


```{r}
mvln = SSdeltaMVLN(base.model1,
                   Fref = "Btgt", 
                   plot = TRUE,
                   addprj=TRUE,
                   virgin = TRUE,
                   mc=100)

```

```{r eval=FALSE}
mvln = SSdeltaMVLN(base.model01,
                   Fref = "Btgt", 
                   plot = TRUE,
                   addprj=TRUE,
                   virgin = TRUE,
                   mc=100)

mvln = SSdeltaMVLN(base.model1,
                   Fref = "Btgt", 
                   plot = TRUE,
                   addprj=TRUE,
                   virgin = TRUE,
                   mc=100)

mvln = SSdeltaMVLN(base.model2,
                   Fref = "Btgt", 
                   plot = TRUE,
                   addprj=TRUE,
                   virgin = TRUE,
                   mc=100)

mvln = SSdeltaMVLN(base.model3,
                   Fref = "Btgt", 
                   plot = TRUE,
                   addprj=TRUE,
                   virgin = TRUE,
                   mc=100)

mvln = SSdeltaMVLN(base.model4,
                   Fref = "Btgt", 
                   plot = TRUE,
                   addprj=TRUE,
                   virgin = TRUE,
                   mc=100)

mvln = SSdeltaMVLN(base.model5,
                   Fref = "Btgt", 
                   plot = TRUE,
                   addprj=TRUE,
                   virgin = TRUE,
                   mc=100)

mvln = SSdeltaMVLN(base.model6,
                   Fref = "Btgt", 
                   plot = TRUE,
                   addprj=TRUE,
                   virgin = TRUE,
                   mc=100)

mvln = SSdeltaMVLN(base.model7,
                   Fref = "Btgt", 
                   plot = TRUE,
                   addprj=TRUE,
                   virgin = TRUE,
                   mc=100)

```

another

```{r eval=FALSE}

nanios<-base.model1$Kobe$Yr
SSBDep<-base.model1$Kobe$B.Bmsy
y_BDPR<-base.model1$Kobe$F.Fmsy


# Integrated Model

kobe1 <- ggplot(base.model1$Kobe)+
  #geom_line(aes(kobebro$BBMSY, kobebro$FFMSY)) +
  geom_rect(aes(xmin = 0, xmax = 0.5, ymin = 0, ymax = 5),
            fill = "#E43338", alpha = 0.5) +
  geom_rect(aes(xmin = 0.5, xmax = 0.75, ymin = 0, ymax = 5),
            fill = "#F2ED23", alpha = 0.5) +
  geom_rect(aes(xmin = 0.75, xmax = 1.25, ymin = 0, ymax = 5),
            fill = "#ACC39A", alpha = 0.5) +
  geom_rect(aes(xmin = 1.25, xmax = 9, ymin = 0, ymax = 1),
            fill = "#608D68", alpha = 0.5) +
  geom_rect(aes(xmin = 1.25, xmax = 9, ymin = 1, ymax = 5),
            fill = "#808080", alpha = 0.5) +
  geom_path(aes(x=base.model1$Kobe$B.Bmsy,y=base.model1$Kobe$F.Fmsy,label=base.model1$Kobe$Yr))+
  geom_point(aes(base.model1$Kobe$B.Bmsy, base.model1$Kobe$F.Fmsy),
             lwd=2) +
  geom_hline(yintercept = 1) +
  geom_vline(xintercept = c(0.5, 0.75, 1.75, 1, 1.25), linetype=2)+
  theme_few()+
  labs(x = expression("BD/BD"[RMS]), y = expression("F/F"[RMS]))+
  # geom_text(data = texto_coords, aes(x = x, y = y, label = etiqueta),
  #            vjust = -0.5)+
  geom_text(aes(x=base.model1$Kobe$B.Bmsy,y=base.model1$Kobe$F.Fmsy,label=base.model1$Kobe$Yr),
            nudge_y = 0.1,size = 3,
            check_overlap = TRUE)
kobe1

```

```{r echo=FALSE}
# model01
tablebias01 <- SShcbias(retroSummary01,quant="SSB",verbose=F)
tablebias01a <- SShcbias(retroSummary01,quant="F",verbose=F)

kbl(tablebias01, booktabs = T,format = "latex",
    caption = "Rho parameter in SSB model s01")  %>% 
    kable_styling(latex_options = "HOLD_position")

kbl(tablebias01a, booktabs = T,format = "latex",
    caption = "Rho parameter in F model s01")  %>% 
    kable_styling(latex_options = "HOLD_position")

# model1
tablebias1 <- SShcbias(retroSummary1,quant="SSB",verbose=F)
tablebias1a <- SShcbias(retroSummary1,quant="F",verbose=F)

kbl(tablebias1, booktabs = T,format = "latex",
    caption = "Rho parameter in SSB model s1")  %>% 
    kable_styling(latex_options = "HOLD_position")

kbl(tablebias1a, booktabs = T,format = "latex",
    caption = "Rho parameter in F model s1")  %>% 
    kable_styling(latex_options = "HOLD_position")

# model2
tablebias2 <- SShcbias(retroSummary2,quant="SSB",verbose=F)
tablebias2a <- SShcbias(retroSummary2,quant="F",verbose=F)

kbl(tablebias2, booktabs = T,format = "latex",
    caption = "Rho parameter in SSB model s2")  %>% 
    kable_styling(latex_options = "HOLD_position")

kbl(tablebias2a, booktabs = T,format = "latex",
    caption = "Rho parameter in F model s2")  %>% 
    kable_styling(latex_options = "HOLD_position")


# model3
tablebias3 <- SShcbias(retroSummary3,quant="SSB",verbose=F)
tablebias3a <- SShcbias(retroSummary3,quant="F",verbose=F)

kbl(tablebias3, booktabs = T,format = "latex",
    caption = "Rho parameter in SSB model s3")  %>% 
    kable_styling(latex_options = "HOLD_position")

kbl(tablebias3a, booktabs = T,format = "latex",
    caption = "Rho parameter in F model s3")  %>% 
    kable_styling(latex_options = "HOLD_position")


# model4
tablebias4 <- SShcbias(retroSummary4,quant="SSB",verbose=F)
tablebias4a <- SShcbias(retroSummary4,quant="F",verbose=F)

kbl(tablebias4, booktabs = T,format = "latex",
    caption = "Rho parameter in SSB model s4")  %>% 
    kable_styling(latex_options = "HOLD_position")

kbl(tablebias4a, booktabs = T,format = "latex",
    caption = "Rho parameter in F model s4")  %>% 
    kable_styling(latex_options = "HOLD_position")


# model5
tablebias5 <- SShcbias(retroSummary5,quant="SSB",verbose=F)
tablebias5a <- SShcbias(retroSummary5,quant="F",verbose=F)

kbl(tablebias5, booktabs = T,format = "latex",
    caption = "Rho parameter in SSB model s5")  %>% 
    kable_styling(latex_options = "HOLD_position")

kbl(tablebias5a, booktabs = T,format = "latex",
    caption = "Rho parameter in F model s5")  %>% 
    kable_styling(latex_options = "HOLD_position")


# model6
tablebias6 <- SShcbias(retroSummary6,quant="SSB",verbose=F)
tablebias6a <- SShcbias(retroSummary6,quant="F",verbose=F)

kbl(tablebias6, booktabs = T,format = "latex",
    caption = "Rho parameter in SSB model s6")  %>% 
    kable_styling(latex_options = "HOLD_position")

kbl(tablebias6a, booktabs = T,format = "latex",
    caption = "Rho parameter in F model s6")  %>% 
    kable_styling(latex_options = "HOLD_position")

# model7
tablebias7 <- SShcbias(retroSummary7,quant="SSB",verbose=F)
tablebias7a <- SShcbias(retroSummary7,quant="F",verbose=F)

kbl(tablebias7, booktabs = T,format = "latex",
    caption = "Rho parameter in SSB model s7")  %>% 
    kable_styling(latex_options = "HOLD_position")

kbl(tablebias7a, booktabs = T,format = "latex",
    caption = "Rho parameter in F model s7")  %>% 
    kable_styling(latex_options = "HOLD_position")

```

## Likelihood

```{r}
like01 <- base.model01$likelihoods_used
like01$model <- rep("s01", nrow(like01))
like01 <- rownames_to_column(like01, var = "Description")
like1 <- base.model1$likelihoods_used
like1$model <- rep("s1", nrow(like1))
like1 <- rownames_to_column(like1, var = "Description")
like2 <- base.model1$likelihoods_used
like2$model <- rep("s2", nrow(like2))
like2 <- rownames_to_column(like2, var = "Description")
like3 <- base.model3$likelihoods_used
like3$model <- rep("s3", nrow(like3))
like3 <- rownames_to_column(like3, var = "Description")
like4 <- base.model4$likelihoods_used
like4$model <- rep("s4", nrow(like4))
like4 <- rownames_to_column(like4, var = "Description")
like5 <- base.model5$likelihoods_used
like5$model <- rep("s5", nrow(like5))
like5 <- rownames_to_column(like5, var = "Description")
like6 <- base.model6$likelihoods_used
like6$model <- rep("s6", nrow(like6))
like6 <- rownames_to_column(like6, var = "Description")
like7 <- base.model7$likelihoods_used
like7$model <- rep("s7", nrow(like7))
like7 <- rownames_to_column(like7, var = "Description")

totalike <- rbind(like01,
                  like1,
                   like2,
                   like3,
                   like4,
                   like5,
                   like6,
                   like7)

```

```{r}
ggplot(totalike)+
  geom_bar( aes(y = reorder(Description, values), x =log(values)),
            stat = "identity", position = "dodge") +
  theme_few() +
  facet_wrap(.~model, 
             ncol=4)+
  labs(y="",
       x="Log Likelihood")+
  theme(axis.text.x = element_text(angle = 90, 
                                   hjust = 1))+
  xlim(0,20)
```

## Likelihood Profile 

```{r eval=FALSE, echo = FALSE}
# 1. Identificar el directorio donde se encuentra el modelo base ----
dirname.model.run <- here("s1")
# 2. Crear un nuevo directorio para el "Perfil_Verosimilitud"  
dirname.R0.profile <- here("Ejercicios","Perfil_Verosimilitud")
dir.create(path=dirname.R0.profile, 
           showWarnings = TRUE, 
           recursive = TRUE)
# 3. Crear un subdirectorio llamado "plots_Verosimilitud" ----
plotdir=paste0(dirname.R0.profile, "/plots_Verosimilitud")
dir.create(path=plotdir,
           showWarnings = TRUE, 
           recursive = TRUE)
# 4. Crear un subdirectorio llamado "simple" ----
reference.dir <- paste0(dirname.R0.profile,'/simple') 
dir.create(path=reference.dir, showWarnings = TRUE, recursive = TRUE)
# 5. Copiar el resultado del modelo base completo en este directorio ----
file.copy(Sys.glob(paste(dirname.model.run, "*.*", sep="/"),
                   dirmark = FALSE),
                    reference.dir)
# 6. Leer la salida del modelo base ----
Base <- SS_output(dir=reference.dir,covar=T)
# 7. Copiar los archivos necesarios de "simple" al directorio "Perfil_Verosimilitud" ----
copy_SS_inputs(dir.old = reference.dir, 
               dir.new =  dirname.R0.profile,
               copy_exe = TRUE,
               verbose = FALSE)
# 8. Leer los archivos del modelo ----
inputs <- r4ss::SS_read(dir = dirname.R0.profile)
# 9. Editar el archivo control la fase de estimación recdev ----
inputs$ctl$recdev_phase <- 1
# 10. Editar el archivo starter para leer los valores de inicio ----
inputs$start$init_values_src <- 0
# 11. Vector de valores para el perfil ----
R0.vec <- seq(15,30,1)  
Nprof.R0 <- length(R0.vec)
# 12. Cambiar el nombre del archivo control en el archivo starter.ss ----
inputs$start$ctlfile <- "control_modified.ss" 
# 13. Incluir prior_like para parámetros no estimados ----
inputs$start$prior_like <- 1                                 
# 14. Escribir los modelos modificados ----
r4ss::SS_write(inputs, dir = dirname.R0.profile, overwrite = TRUE)
# 15. Ejecutar la función profile() ----
#?SS_profile()
profile <- profile(dir=dirname.R0.profile, # directory
                      exe="ss_osx",
                      oldctlfile ="control.ss",
                      newctlfile="control_modified.ss",
                      string="SR_LN(R0)",
                      profilevec=R0.vec)
# 16. Leer los archivos de salida ----
# (con nombres como Report1.sso, Report2.sso, etc.)
prof.R0.models <- SSgetoutput(dirvec=dirname.R0.profile, 
                              keyvec=1:Nprof.R0, 
                              getcovar = FALSE) 
# 17. Resumir las salidas con la función SSsummarize()  ----
prof.R0.summary <- SSsummarize(prof.R0.models)
# 18. Identificar los componentes de Verosimilitud ----
mainlike_components         <- c('TOTAL',
                                 "Survey", 
                                 'Length_comp',
                                 "Age_comp",
                                 "Catch",
                                 'Size_at_age',
                                 'Recruitment') 
mainlike_components_labels  <- c('Total likelihood',
                                 'Index likelihood',
                                 'Length likelihood',
                                 "Age likelihood",
                                 "Catch Likelihood",
                                 'Size_at_age likelihood',
                                 'Recruitment likelihood') 

```

Funciones para generar plots de perfil de verosimilitud ----

### SSplotProfile() ----

```{r eval=FALSE, echo = FALSE}
png(file.path(plotdir,"R0_profile_plot.png"),
    width=7,
    height=4.5,
    res=300,
    units='in')
par(mar=c(5,4,1,1))
SSplotProfile(prof.R0.summary,           # summary object
              profile.string = "R0",     # substring of profile parameter
              profile.label=expression(log(italic(R)[0])), 
              ymax=2050,minfraction = 0.001,
              pheight=4.5, 
              print=FALSE, 
              plotdir=plotdir, 
              components = mainlike_components, 
              component.labels = mainlike_components_labels,
              add_cutoff = TRUE,
              cutoff_prob = 0.95)

Baseval <- round(Base$parameters$Value[grep("R0",Base$parameters$Label)],2)
Baselab <- paste(Baseval,sep="")
axis(1,at=Baseval,label=Baselab)
abline(v = Baseval, lty=2)
dev.off()
```

## Convergence criteria

0.0001 final convergence criteria (e.g. 1.0e-04) 


```{r eval=FALSE, echo = FALSE}
# Comparación de series de tiempo 
labs <- paste("SR_Ln(R0) = ",R0.vec)
labs[which(round(R0.vec,2)==Baseval)] <- paste("SR_Ln(R0) = ",
                                               Baseval,"(Base model)")

SSplotComparisons(prof.R0.summary,
                  legendlabels=labs,
                  pheight=4.5,png=TRUE,
                  plotdir=plotdir,
                  legendloc='bottomleft')


```

piner Plot

```{r echo = FALSE}
#### R0_profile_plot_Length_like ----
png(file.path(plotdir,"R0_profile_plot_Length_like.png"),
    width=7,
    height=4.5,
    res=300,
    units='in')
par(mar=c(5,4,1,1))
PinerPlot(prof.R0.summary, 
          profile.string = "R0", 
          component = "Length_like",
          main = "Changes in length-composition likelihoods by fleet",
          add_cutoff = TRUE,
          cutoff_prob = 0.95)
Baseval <- round(Base$parameters$Value[grep("SR_LN",
                                      Base$parameters$Label)],2)
Baselab <- paste(Baseval,sep="")
axis(1,at=Baseval,
     label=Baselab)
abline(v = Baseval, lty=2)
dev.off()
```


```{r echo = FALSE}
#### R0_profile_plot_Survey_like ----
png(file.path(plotdir,"R0_profile_plot_Survey_like.png"),
    width=7,
    height=4.5,
    res=300,
    units='in')
par(mar=c(5,4,1,1))
PinerPlot(prof.R0.summary, 
          profile.string = "R0", 
          component = "Surv_like",
          main = "Changes in Index likelihoods by fleet",
          add_cutoff = TRUE,
          cutoff_prob = 0.95, legendloc="topleft")
Baseval <- round(Base$parameters$Value[grep("SR_LN",
                                            Base$parameters$Label)],2)
Baselab <- paste(Baseval,sep="")
axis(1,at=Baseval,label=Baselab)
abline(v = Baseval, lty=2)
dev.off()
```

## Outputs Model Base

```{r}
outps1 <- base.model1$timeseries[1:53, 2:8]

# addWorksheet(wb, "variable_s2")
# writeData(wb, "variable_s2", outps1)
# 
# # Guardar el workbook
# saveWorkbook(wb, "DataKrill.xlsx", overwrite = TRUE)
# 
# 
# out_s2 <- read_excel("DataKrill.xlsx", 
#     sheet = "variable_s2")
outps1 %>%
  kbl(booktabs = TRUE,
      format = "html",
    caption = "Main variables outputs from stock asssessment krill in WAP") %>%
  kable_paper("hover", 
              full_width = TRUE)%>%
  kable_styling(latex_options = c("striped",
                                  "condensed"),
                full_width = FALSE,
                font_size=9,
                html_font = "arial")#%>% 
  #pack_rows(index = c("Estimation" = 1,
   #                     "Prediction" = 45))
```

## Comparision outputs betwwen scenarios

```{r echo=FALSE}

#read in all model runs
#note if cover=T you need a hessian; if covar=F you do not need a hessian
biglist <- SSgetoutput(keyvec = NULL,
                       dirvec = c(dir01,
                                  dir1,
                                  dir2,
                                  dir3,
                                  dir4,
                                  dir5,
                                  dir6,
                                  dir7),
                       getcovar = F)

#create summary of model runs from list above
summaryoutputall <- SSsummarize(biglist)
```


```{r echo=FALSE}
compall <- SSplotComparisons(summaryoutputall,
                  legendlabels = c("No Spatial",
                                   "s1",
                                   "s2",
                 "s3",
                 "s4",
                 "s5",
                 "s6",
                 "s7"),
                 labels = c("Year", 
                            "Spawning biomass (t)",
                            "Relative spawning biomass", 
                            "Age-0 recruits (1,000s)",
                            "Recruitment deviations", 
                            "Index", "Log index", 
                            "1 - SPR", 
                            "Density",
                            "Management target", 
                            "Minimum stock size threshold",
                            "Spawning output",
                            "Harvest rate"),
                 png=TRUE,plotdir=Figs)

```

```{r}
comtable <- SStableComparisons(summaryoutputall,
                   likenames = c("TOTAL", 
                                 "Survey", 
                                 "Length_comp",
                                 "Age_comp", 
                                 "priors",
                                 "Size_at_age"), 
                   names = c("Recr_Virgin",
                             "R0", 
                             "steep",
                             "NatM",
                             "L_at_Amax", 
                             "VonBert_K", 
                             "SSB_Virg", 
                             "Bratio_2020",
                             "SPRratio_2020"),
                   digits = NULL,
                   modelnames = c("s01", "s1","s2",
                 "s3",
                 "s4",
                 "s5",
                 "s6",
                 "s7"),
                   csv = F,
                   csvdir =~"/DOCAS/SA_Krill",
                   csvfile = "parameter_comparison_table.csv",
                 verbose = TRUE,
                   mcmc = FALSE)
kbl(comtable, booktabs = T,format = "latex",
    caption = "Comparacion likelihood y parámetros s01, s1, s2, s3, s4, s5, s6 y s7")  %>% 
    kable_styling(latex_options = "scale_down")

```

## Comparision between select models `No spatial`and `Spatial implicit` and `Spatial W/ new set parametres`

```{r echo=FALSE}
#read in all model runs
#note if cover=T you need a hessian; if covar=F you do not need a hessian
biglist01_1_7 <- SSgetoutput(keyvec = NULL,
                       dirvec = c(
                                  dir01,
                                  dir1,
                                  dir7),
                       getcovar = F)

#create summary of model runs from list above
summaryoutput01_1_7 <- SSsummarize(biglist01_1_7)

SSplotComparisons(summaryoutput01_1_7,
                  legendlabels = c(
                                   "No Spatial",
                 
                 "Spatial",
                 "Spatial w/ new set LH paramters"),
                 filenameprefix = "COM1",
                 labels = c("Year", 
                            "Spawning biomass (t)",
                            "Relative spawning biomass", 
                            "Age-0 recruits (1,000s)",
                            "Recruitment deviations", 
                            "Index", "Log index", 
                            "1 - SPR", 
                            "Density",
                            "Management target", 
                            "Minimum stock size threshold",
                            "Spawning output",
                            "Harvest rate"),
                 png=TRUE,plotdir=Figs)
```

comparision between select models `Old Paramters` and `New Parameters WG SAM 2024/23`

```{r echo=FALSE}
#read in all model runs
#note if cover=T you need a hessian; if covar=F you do not need a hessian
biglist1_7 <- SSgetoutput(keyvec = NULL,
                       dirvec = c(
                                  dir1,
                                  dir7),
                       getcovar = F)

#create summary of model runs from list above
summaryoutput1_7 <- SSsummarize(biglist1_7)

comp1_7 <- SSplotComparisons(summaryoutput1_7,
                  legendlabels = c(
                                   "Old Parameters",
                 
                 "Parameters WG SAM 2024/27"),
                 labels = c("Year", 
                            "Spawning biomass (t)",
                            "Relative spawning biomass", 
                            "Age-0 recruits (1,000s)",
                            "Recruitment deviations", 
                            "Index", "Log index", 
                            "1 - SPR", 
                            "Density",
                            "Management target", 
                            "Minimum stock size threshold",
                            "Spawning output",
                            "Harvest rate"),
                 png=TRUE,plotdir=Figs)
```
### Comparision between select models `s1` and `s9`

```{r echo=FALSE}
#read in all model runs
#note if cover=T you need a hessian; if covar=F you do not need a hessian
biglist1_9 <- SSgetoutput(keyvec = NULL,
                       dirvec = c(
                                  dir1,
                                  dir9),
                       getcovar = F)

#create summary of model runs from list above
summaryoutput1_9 <- SSsummarize(biglist1_9)

comp1_9 <- SSplotComparisons(summaryoutput1_9,
                  legendlabels = c(
                                   "Base",
                 "selectivity Block 2006"),
                 filenameprefix = "COM1_9",
                 labels = c("Year", 
                            "Spawning biomass (t)",
                            "Relative spawning biomass", 
                            "Age-0 recruits (1,000s)",
                            "Recruitment deviations", 
                            "Index", "Log index", 
                            "1 - SPR", 
                            "Density",
                            "Management target", 
                            "Minimum stock size threshold",
                            "Spawning output",
                            "Harvest rate"),
                 png=TRUE,plotdir=Figs)
```
## Comparsion in sd long term time series

```{r fig.height=6, fig.width=8}

sblt01 <- base.model01$timeseries$Bio_all[1:50]
sblt1 <-  base.model1$timeseries$Bio_all
sblt2 <-  base.model2$timeseries$Bio_all
sblt3 <-  base.model3$timeseries$Bio_all
sblt4 <-  base.model4$timeseries$Bio_all
sblt5 <-  base.model5$timeseries$Bio_all
sblt6 <-  base.model6$timeseries$Bio_all
sblt7 <-  base.model7$timeseries$Bio_all

errt <- as.data.frame(cbind(sblt01,
              sblt1,
              sblt2,
              sblt3,
              sblt4,
              sblt5,
              sblt6,
              sblt7))
err_long <- errt %>%
  pivot_longer(cols = everything(),
               names_to = "variable",
               values_to = "value")  
titles <- c("No Spatial",
            "Spatial implicit", 
            "Spatial w/ without S-R relation",
            "Spatial w/ S-R relation weak (0.9 steepness)",
            "Spatial w/ S-R relation strong (0.6 steepness)",
            "Spatial w/ S-R relation mid-strong estimated",
            "Spatial w/ Ricker S-R relation estimated", 
            "Spatial w/ New set parameters Linf, k and M")


ggplot(err_long, aes(x = variable, y = value)) +
  geom_boxplot() +
  geom_jitter(width = 0.2, alpha = 0.3) +
  labs(title = "Long Term Predictions",
       x = "",
       y = "Spawning Biomass (t)") +
  theme_few() +
  scale_x_discrete(labels= titles)+
  theme(axis.text.x = element_text(angle = 90,
                                   hjust = 1))
```

## Autocorrelation in Recruit and Biomas

```{r fig.width=8, fig.height=7}
recs01 <- base.model01$recruit$pred_recr
recs1 <- base.model1$recruit$pred_recr
recs2 <- base.model2$recruit$pred_recr
recs3 <- base.model3$recruit$pred_recr
recs4 <- base.model4$recruit$pred_recr
recs5 <- base.model5$recruit$pred_recr
recs6 <- base.model6$recruit$pred_recr
recs7 <- base.model7$recruit$pred_recr

p01 <- ggAcf(base.model01$recruit$pred_recr) +
  ggtitle("No Spatial") +
  theme_few()
p1 <- ggAcf(base.model1$recruit$pred_recr) +
  ggtitle("Spatial implicit") +
  theme_few()
p2 <- ggAcf(base.model2$recruit$pred_recr) +
  ggtitle("Spatial w/ without S-R relation") +
  theme_few()
p3 <- ggAcf(base.model3$recruit$pred_recr) +
  ggtitle("Spatial w/ S-R relation weak (0.9 steepness)") +
  theme_few()
p4 <- ggAcf(base.model4$recruit$pred_recr) +
  ggtitle("Spatial w/ S-R relation strong (0.6 steepness)") +
  theme_few()
p5 <- ggAcf(base.model5$recruit$pred_recr) +
  ggtitle("Spatial w/ S-R relation mid-strong estimated") +
  theme_few()
p6 <- ggAcf(base.model6$recruit$pred_recr) +
  ggtitle("Spatial w/ Ricker S-R relation estimated") +
  theme_few()
p7 <- ggAcf(base.model7$recruit$pred_recr) +
  ggtitle("Spatial w/ New set parameters Linf, k and M") +
  theme_few()




ggarrange(p01,p1,
          p2,
          p3,
          p4,
          p5,
          p6,
          p7, ncol=3, nrow = 3,
          common.legend = T)

```

## Recruit deviation

```{r}
dev01 <- base.model01$recruit[1:43,c(1,7)] %>% 
  mutate(Serie = "S01")
dev1 <- base.model1$recruit[1:43,c(1,7)]%>% 
  mutate(Serie = "S1")
dev2 <- base.model2$recruit[1:43,c(1,7)] %>% 
  mutate(Serie = "S2")
dev3 <- base.model3$recruit[1:43,c(1,7)] %>% 
  mutate(Serie = "S3")
dev4 <- base.model4$recruit[1:43,c(1,7)] %>% 
  mutate(Serie = "S4")
dev5 <- base.model5$recruit[1:43,c(1,7)] %>% 
  mutate(Serie = "S5")
dev6 <- base.model6$recruit[1:43,c(1,7)] %>% 
  mutate(Serie = "S6")
dev7 <- base.model7$recruit[1:43,c(1,7)] %>% 
  mutate(Serie = "S7")

data_list <- list(dev01, dev1, dev2, dev3, dev4, dev5, dev6, dev7)
titles <- c("No Spatial",
            "Spatial implicit", 
            "Spatial w/ without S-R relation",
            "Spatial w/ S-R relation weak (0.9 steepness)",
            "Spatial w/ S-R relation strong (0.6 steepness)",
            "Spatial w/ S-R relation mid-strong estimated",
            "Spatial w/ Ricker S-R relation estimated", 
            "Spatial w/ New set parameters Linf, k and M")
# Lista para almacenar los gráficos
plots <- list()

# Crear los gráficos en un ciclo for
for (i in 1:length(data_list)) {
  p <- ggplot(data_list[[i]], aes(x = Yr, y = dev)) +
    geom_point() +
    stat_smooth(method = "loess", 
                span = 0.15, 
                se = T) +
    geom_hline(yintercept = 0,
               col="red",
               linetype = "dotdash")+
    ggtitle(titles[i]) +
    theme(axis.text.x = element_text(angle = 90,
                                   hjust = 1))+
    theme_few()
  plots[[i]] <- p
}

# Organizar los gráficos con ggarrange
combined_plot <- ggarrange(plotlist = plots, ncol = 3, nrow = 3)


combined_plot
```
## Platoons analisis

```{r}
# Definir los valores del eje x
x_values <- seq(0, 7, by=0.1)

# Definir las medias y desviaciones estándar para los valores específicos
means <- c( 2.9,  3.5, 4.1)
std_devs <- c(0.65, 0.8, 0.65)


# Crear un dataframe que contenga los valores de x, la media y las desviaciones estándar
data <- data.frame(
  x = rep(x_values, each = length(means)),
  mean = rep(means, times = length(x_values)),
  std_dev = rep(std_devs, times = length(x_values))
)

# Calcular las curvas de las desviaciones estándar
data <- data %>%
  mutate(y = exp(-((x - mean)^2) / (2 * std_dev^2)))

# Graficar con ggplot
ggplot(data, aes(x = x, y = y, linetype = as.factor(mean))) +
  geom_line(size = 1) +
  scale_linetype_manual(values = means,
                        name ="Platoon") +
  labs(x = "", y = "growth increment") +
  theme_minimal()+
  xlim(0,7)

```
## AKL

In a catch-at-length model like krill assessment the AKL matrix is modelled trought parametrization process

```{r AKL, out.width='60%',  fig.cap="Representation of ALK Matrix to krill in 48.1"}
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 = "True Age",
       y = "Length",
       fill = "Value") +
  theme_bw()+
  theme(legend.position = "none")
```


\newpage

# REFERENCES